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Abstract 

We give an overview of the fundamental concepts of density functional 
theory. We give a careful discussion of the several density functionals and 
their differentiability properties. We show that for nondcgcncratc ground 
states we can calculate the necessary functional derivatives by means of 
linear response theory, but that there arc some differentiability problems 
for degenerate ground states. These problems can be overcome by extending 
the domains of the functionals. We further show that for every interacting 
w-representablc density we can find a noninteracting u-rcprcscntable density 
arbitrarily close and show that this is sufficient to set up a Kohn-Sham 
scheme. We finally describe two systematic approaches for the construction 
of density functionals. 
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1. Introduction 



In this paper we give an overview of the foundations of density functional 
theory for stationary systems. In the discussion we try to be as precise 
as possible and we therefore carefully discuss several exact density func- 
tional and the domains of densities on which they are defined. At the 
heart of almost any application of density functional theory lie the Kohn- 
Sham equations. These equations describe a noninteracting system that 
is required to produce the ground state density of an interacting system. 
Therefore any discussion of the validity of these equations has to focus 
on how well we can approximate a density of an interacting system with 
that of a noninteracting system. Closely related to this question is whether 
or not the exact functionals that we defined have functional derivatives. 
This is because the various potentials, such as the Kohn-Sham potential, 
in density functional theory appear as functional derivatives of the energy 
with respect to the density. The investigation of the existence of functional 
derivatives will form a central theme of this review. 

A large part of this work will follow the proofs of the fundamental papers 
by Licb [1] and Englisch and Englisch [2, 3]. In this paper we try to make 
these two important works more accessible by providing some background 
on the mathematics involved. We further show, from a more physical view- 
point, how to calculate functional derivatives by response theory and show 
that for nondegenerate ground states the static density response function is 
invertable. We also show that for every interacting w-representable density 
there is a noninteracting w-representable density arbitrarily close to it and 
that this is sufficient to set up a Kohn-Sham scheme. We finally discuss 
two systematic approaches for the construction of the exchange-correlation 
functional and provide an outlook and conclusions. 

2. Definition of the problem 

Density-functional theory in its earliest formulation by Hohenbcrg, Kohn 
and Sham [5, 8] aims at a description of the description of ground state 
properties of many-electron systems in terms of the electron density. One 
may wonder why this is possible. Let us therefore investigate this point 
more closely. Consider a Hamiltonian of a stationary many-body system 

H v = f+V + W (1) 

where T is the kinetic energy of the electrons, V the external potential, 
and W the two-particle interaction. We denote the Hamilton operator H v 
with a subindex v to indicate that we will consider the Hamiltonian as a 
functional of the external potential v. The constituent terms are explicitly 



given as 



T = 



N 1 

( 2 ) 



2 

i=i 

N 



V = 2«( r< ) (3) 

i=\ 
N 

W = ^mr.-r.l) (4) 

i>j 

where w(|r|) = l/|r| will in our discussion always be the Coulomb potential. 
We are interested in electronic systems, i.e. molecules and solids. For all 
these systems the kinetic energy operator T and two-particle interaction 
W are identical. They only differ in the form of the external potential v(r) 
and the number of electrons N. The properties of all these systems can 
therefore be regarded as a functional of the external potential v. This is in 
particular the case for the ground state wavefunction ^[t 1 ]) and the ground 
state energy E[v], which are of course related by the Schrodinger equation 

(f+V + W)\*[v])=E[v]\*[v]) (5) 

At first sight we have not gained much by this viewpoint, the problem looks 
as difficult as before. However, the problem will have a different appearance 
once we eliminate the potential in favor of the electron density. One may 
wonder what is so particular about the relation between the density and 
the potential. Let us therefore look specifically at the term which describes 
the external potential. It can be written as 



V = / cfYw(r)n(r) (6) 



where we defined the density operator by 

JV 



i(r) = 2«y(r- ri ) (7) 



»=i 

The expectation value of V is given by 



(tflVI*) = J d 3 rn(r)v(r) (8) 



where n(r) the electron density . The electron density is obtained from 
the many-body wavefunction, which in our case will usually be the ground 



state wavefunction of Hamiltonian H v , by 



n(r 1 ) = <*|n(r 1 )|*)=JV ^ / d 3 r 2 . . . d 3 r N |f (naj, . . . , r N a N )\ 2 (9) 



In this expression is the spin variable for electron i. The physical inter- 
pretation of the density is that n(r)AV is proportional to the probability 
of finding an electron in an infinitesimal volume AV around r. Note that 
in the definition Eq.(9) we sum over the spin coordinate o\ so we do not 
distinguish between the probabilities of finding up or down spin electrons 
at r. We see from Eq.(8) that the external potential v(r) and the electron 
density n(r) are conjugate variables. With this we mean that they occur 
as a simple product in the contribution of the external potential to the 
energy. It is exactly this property that allows us to prove the Hohenberg- 
Kohn theorem which establishes a 1-1-correspondence between the density 
and the external potential. We can therefore go from a functional E[v] of 
the external potential to a functional E[n] of the density. We can ask our- 
selves whether there is any reason that E[n] would be easier to calculate 
than E[v\. A priori there is no reason to expect this. However, we will 
see that the determination of E[n] is equivalent to the solution of a set of 
one-particle equations, known as the Kohn-Sham equations, with a poten- 
tial v s [n] that is also a functional of the density. By now we know that 
we can find practical and useful approximations for this potential v s [n] and 
that the Kohn-Sham equations have been succesfully applied to the calcu- 
lation of properties of many electronic systems. However, as in any physical 
theory, there are a number of assumptions made in the transformation to 
the Kohn-Sham one-particle equations. The aim of this paper is to discuss 
these assumptions and thereby the validity of the Kohn-Sham equations. 

3. Conditions on the electron density and the external 
potentials 

We begin by giving a discussion of the properties of the two key objects in 
density functional theory, the density n and the external potential v. The 
density has the obvious properties 



These properties follow directly from the definition of the density and the 
usual normalization condition on the wavefunction. If we take into account 
that the density is obtained as the density of a bound eigenstate of Hamil- 
tonian (1) we can derive further conditions. For this we put the physical 



<7l...<7jV 



n(r) > 




(10) 



constraint on the many-body system that it has a finite expectation value 
of the kinetic energy, i.e. 



N 

T M=2E E J d 3 r 1 ...d 3 r N \V i *\ 2 <™ (11) 



i — 1 (7 1 ...fT jv 



At this point it is useful to introduce a new space of functions. We say that 
a function / is in H 1 (TZ n ) (1Z denotes the real numbers) if 

d"r(l/W| 2 + |V/(r)| 2 y /2 <^ (12) 

The space of functions H 1 (TZ n ) is called a Sobolev space. The supindex 1 
refers to the fact that the definition of the norm contains only first order 
derivatives. We therefore see that finiteness of the kinetic energy implies 
that is an element of the function space i/ 1 (7?. 3Ar ). Differentiation of 
Eq.(9) and use of the Cauchy-Schwarz inequality then leads to [1] 

(Vm(n)) 2 < 4JVn(n) J2 J ^ ■ ■ ■ rf 3 rjv|Vi*| 2 (13) 

and hence 

l - j d 3 r(V\AiW) 2 < T[tf] < oo (14) 

We see that the finiteness of the kinetic energy puts a constraint on the 
density. From Eq.(12) we see that y/n belongs to i? 1 (7?. 3 ). We also see that 
if we consider systems with a finite kinetic energy, then we only need to 
consider the following set of densities 

S = {n\ n(r) > 0, \fn G H 1 ^ 3 ), J d 3 rn(r) = N} (15) 

This set of densities has a property which will be of importance later, 
namely S is convex. With this we mean that if ri\ and n 2 arc elements of 
S, then also Xn\ + (1 — A)n 2 is an clement of S where < A < 1. This 
property is easily proven using the Cauchy-Schwarz inequality. Now we will 
derive constraints on the allowed set of external potentials. In order to do 
this we introduce some other function spaces. We say that a given function 
/ belongs to the space L p if 



(/ 



i/p 

d 3 r\f(v)\P) <^ (16) 



Note that here we only consider functions on the usual thrccdimcnsional 
coordinate space 1Z 3 . The letter L refers to Lebesque integration, a feature 



that assures that the function spaces are complete (complete normed spaces 
are also called Banach spaces). We will, however, not go into the detailed 
mathematics and refer the interested reader to the literature [4]. We just 
note that for continuous functions the integral is equivalent to the usual 
(Ricmann) integral. Equation (16) defines a norm on the space L p and 
we see from Eq.(lO) that the density belongs to L 1 . From the condition 
of finite kinetic energy and the use of a Sobolev inequality one can show 
that [1] 

J d 3 rn 3 (r) <C J d 3 r(vVn(r)) 2 < oo (17) 

where C — 3(7r/2) 4 / 3 . In other words, the finiteness of the kinetic energy 
implies that the density is also in the space L 3 . Since we already know 
that n <G L 1 we find that the density is clement of the intersection of both 
spaces, i.e. n € L 1 n L 3 . We also see from the inequality (17) that S is a 
subset of L 1 n L 3 . Let us now see what this implies for the allowed set of 
external potentials when we require the expectation value of the external 
potential to be finite, i.e. we require 



d 3 rn(r)v(r) 



<||nu||i<oo (18) 



If the potential is bounded, i.e. |w(r)| < M for some finite number M then 



/ 



d 3 rn(r)v(r) 



< sup |w(r)|iV < oo (19) 



The space of bounded functions is called L°° and has the norm 

H/Hoo =supj/(r)| (20) 

The supremum is defined to be the smallest number M such that |/(r) | < M 
almost everywhere. The term 'almost everywhere' has a precise mathemat- 
ical meaning for which we refer to the literature [4]. We almost never use 
it in the remainder of this paper. We therefore conclude that if v € L°° 
then the expectation value of the external potential is finite. To show this 
we used that n £ L 1 . But we also know that n € L 3 and if we make use of 
the Holder inequality 

\\.fgh<\\f\\ P \\g\\ q (21) 

with 1/p+l/q = 1 we obtain 

H™||i < IMUIMIf (22) 

which is finite if v € L § . Therefore the most general set of potentials for 
which the expectation value (^IV^I^) is finite, is the set 

+L X = {v\v = u + w,ue L^,w e L°°} (23) 



i.e. the set of potentials that can be written as a sum of a function from L 2 
and a function from L°° . This is also a normed function space with norm 

||v||=inf{||u||| + ||u'|| 00 |w = it + w} (24) 

In the remainder of this paper we will always consider the densities to be in 
the space L 1 n L 3 and the potentials in the space L 't- + L°° . It is important 
to note that the Coulomb potential is in the latter set since we can write 

l = g(i-H) + g(H-p (25) 

|r| |r| |r| 

where is the Heaviside function, 0{x) = if x < and 0{x) = 1 if 
x > 0. One can readily check that the first and the second term on the 
right hand side are in is and L°° respectively. One can easily extend this 
result to a finite sum of Coulomb potentials and therefore every molecule 
can be described with the potentials in the space L 2 + L°° . One may finally 
wonder what the condition of finite Coulombic electron-electron repulsion, 
i.e. (\E'|W|\I/) < 00, would imply for the density. However, one can show 
that if the kinetic energy (*|T|*) is finite then also (*|W'|*) is finite [1], 
so this does not yield any new constraints on the density. 



4. The Hohenberg-Kohn theorem 

The basis of density-functional theory is provided by the Hohenberg-Kohn 
theorem [5]. We will provide a proof of this theorem for nondegenerate 
ground states. The case of degenerate ground states will be discussed later. 
The Hohenberg-Kohn theorem states that the density n(r) of a nondegen- 
erate ground state uniquely determines the external potential v(r) up to an 
arbitrary constant. This means that the external potential is a well-defined 
functional w[n](r) of the density. 

In the proof essential use is made of the fact that the density and the po- 
tential are conjugate variables. For the same reason we can, for instance, 
prove that the 2-particle interaction is a unique functional of the diago- 
nal 2-particle density matrix. The general mapping between ./V-particle 
density matrices and TV-body potentials is discussed by De Dominicis and 
Martin [6]. 

3 

Let us consider the subset VdH L°° of potentials that yield a normaliz- 
able nondegenerate ground state. The solution of the Schrodingcr equation 
provides us with a mapping from the external potential to the ground state 
wavefunction, v(r) — > ^[w])- Since we assume that we are dealing with 
nondegenerate ground states is uniquely determined apart from a 

trivial phase factor. We have therefore established a map C : V — > <&, 
where <J> is the set of ground states. 



We will first prove that the map C is invertible. Suppose that and 
1^2) G _ff 1 (7?. 3Ar ) correspond to external potentials v\ and v 2 G L°° + 
where V\ 7^ v 2 + C. We have to show that 7^ l 1 !^)- If we assume that 
l^i) = 1*2) = I*) then by subtraction of the Hamiltonian for and 
1*2) we find that 

(Vi-Va)!*) - (Ei -E 2 )\V) (26) 

If v\ — v 2 is not constant in some region then \l/ must vanish in this region 
for the above equation to be true. However if v\,v 2 G L°° + then |^) 
cannot vanish on an open set (a set with nonzero measure) by the unique 
continuation theorem [1]. So we obtain a contradiction and hence we must 
have made a wrong assumption. Therefore 7^ |W 2 ) and we obtain the 
result that different potentials ( differing more than a constant) give differ- 
ent wavefunctions. Consequently we find that the map C is invertible. 
We now define the set A as the set of densities which come from a nonde- 
generate ground state, where we only consider ground state densities from 
potentials in the set L°° + . The set A is obviously a subset of the previ- 
ously defined set S. From a given wavefunction in the set of ground states 
<I> we can calculate the density according to Eq.(9). This provides us with 
a second map D : Q ^ A from ground state wavefunctions to ground state 
densities. Also this map is invertible. To show this we calculate 

E[ Vl ] = (*H|f +Vi + W1*M) < <*H|f +Vi + w\y[v 2 }) 
= (*H|f +V 2 + W]*H) + J n 2 (r)( Vl (r) - v 2 (r))dr 

= E[v 2 } + Jn 2 (r)( Vl (r)-v 2 (r))dr (27) 

Likewise we find 

E[v 2 ] < E[ Vl ] + J m(r)(«2(r) - «i(r))dr (28) 

Adding both inequalities then yields the new inequality 

d 3 r(n 2 (r) - ni (r))(v 2 (r) - Vi (t)) < (29) 

If we assume that m — n 2 then we obtain the contradiction < and we 
conclude that different ground states must yield different densities. There- 
fore the map D is also invertible. Consequently the map DC : V — > A is 
also invertible and the density therefore uniquely determines the external 
potential. This proves the Hohcnbcrg-Kohn theorem. 

Let us now pick an arbitrary density out of the set A of densities of non- 
degenerate ground states. The Hohenberg-Kohn theorem then tells us that 



there is a unique external potential v (to within a constant) and a unique 
ground state wavefunction (to within a phase factor) corresponding 

to this density. This also means that the ground state expectation value of 
any observable, represented by an operator O, can be regarded as a density 
functional 

0[n] = (tf[n]|6|tf[n]) (30) 

In particular we can thus define the Hohcnbcrg-Kohn functional Fhk on 
the set A as 

F HK [n] = ^>[n\\f + W\^[n\) (31) 
With this functional we can define the energy functional E v as 

E v [n] = J n{r)v{r)dr + F HK [n] (32) 

If n is a ground state density corresponding to external potential v and 
n an arbitrary other ground state density then 

E VQ [n] = Jn{r)v {r)dr + F HK [n] = (^[n]\f+% + W\^[n]) 

> (*H|f + V + Wmn ])=E Vo [n ] (33) 

Therefore 

E[v] = in^jy n(r)v(r)dr + F HK [n]^ (34) 

and we see that the ground state energy of a many-body system can be 
obtained by minimization of a density functional. For application of this 
formula we have to know Fhk on the set A. In practice it is, of course, 
impossible to calculate Fhk exactly on this set of densities. Instead one 
would prefer to make an explicit approximation for Fhk as is usually done 
within the Kohn-Sham scheme. Before we go into that let us first discuss 
some properties of Fhk- The functional Fhk is a convex functional, i.e. if 
ni,n 2 € A and if Aini + \2ri2 € A with < Ai, A 2 < 1 and Ai + A 2 = 1 
then 

F HK [Mni + X 2 n 2 ] < \xFhk\ui] + X 2 F H K[n 2 ] (35) 

This is readily proven. Suppose that the ground state densities ni, n 2 , Ai7ii + 
A 2 n 2 e A correspond to the external potentials v\,v 2 and v. Then 

F HK [n] = (V[n\\f + V + W\V[n\)- j d 3 rn(r)v(r) 

= Ai<*[n]|fl"|*[n]) + A 2 <*[n]|ff|*[n]) - J d 3 rn(r)v(r) 
< Ai<*[m]|f + Wmm]) + A 2 (*[n 2 ]|f + W\^[n 2 \) 



+ J (Aini(r) + A 2 n 2 (r))i;(r)cf 3 r - J d 3 rn(r)v(r) 

= XiF HK [n 1 } + X 2 F HK [n 2 } (36) 

and we obtain the convexity of Fhk- Note, however, that the domain A 
of Fhk docs not need to be convex, i.e. if n\,ri2 € A then not necessarily 
AiTii + A 2 n 2 € A with < Ai, A 2 < 1 and \i + A 2 = 1. We will come back 
to this point later when we consider the differentiability of Fhk- Let us 
first collect our results in the form a theorem 

Theorem 1 (Hohenberg-Kohn) The density n corresponding to a non- 
degenerate ground state specifies the external potential v up to a constant 
and the ground state wavefunction \^[n]) up to a phase factor. Moreover, 

1. Any ground state expectation value corresponding to an observable O 
is a functional of the density according to 

0[n] = (*[n]|6|*[n]) 

2. The ground state energy of a system with a nondegenerate ground 
state and an external potential v can be obtained from 



E[v}= inf <^ / n(r)v(r)dr + F HK \n] 

neA I J 

where F HK [n} = (V[n]\f + W\V[n]) . 
3- F H k is convex 

We finally discuss an interesting consequence of the Hohenberg-Kohn the- 
orem. Suppose that in Eq.(29) we take v 2 = v± + eSv where 5v is not a 
constant. We therefore obtain potential w 2 by a small perturbation from 
potential vi. By means of perturbation theory we can then calculate the 
change in the density which gives 

n 2 (r) =n 1 (r) + e<5n(r) + 0(e 2 ) (37) 

where Sn can be calculated from the static density response function \ 

Sn(r) = J d 3 r' X (r,r')Sv(r') (38) 

The properties and explicit form of x are described in more detail in a later 
section. Therefore 

Mr) - ^(r))(n 2 (r) - m(r)) = e 2 8v{v)5n{v) + 0(e 3 ) (39) 



If we insert this expression in Eq.(29) and divide by e 2 > then we obtain: 



J d 3 rSv(r)6n(r) + 0(e) < (40) 

Now taking the limit e — > and expressing 8n in terms of the response 
function we obtain 

J d 3 rd 3 r'6v{r)x{r,r')5v{r') < (41) 

This is true for an arbitrary nonconstant potential variation. We therefore 
see that the eigenvalues of x, when we regard x as an integral operator, 
are negative. Moreover we see that the only potential variation that yields 
a zero density variation is given by Sv = C where C is a constant. This 
implies that x 1S invcrtablc. We will go more closely into this matter in a 
later section where we will prove the same using the explicit form of the 
density response function. 



5. Kohn-Sham theory by Legendre transforms 

The method described in this section goes back to the work of De Do- 
minicis and Martin [6] . This work discusses the relations between iV-body 
potentials and iV-particle density matrices, of which the density-potential 
relation to be discussed here is a special case. The ground state energy E[v] 
and wave function \^f[v]) are considered to be functionals of the external 
potential through solving the time-independent Schrodinger equation 

(f+V + W)\^[v])=E[v]\^[v\) (42) 

where the two-particle interaction W is kept fixed. From this equation we 
see that the ground state energy as a functional of the external potential v 
can also be written as 

E[v] = (9[v]\H v \*[v]) (43) 

Our goal is now to go from the potential as our basic variable, to a new 
variable, which will be the electron density. The deeper reason that this is 
possible is that the density and the potential are conjugate variables. With 
this we mean that the contribution of the external potential to the total 
energy is simply an integral of the potential times the density. We make use 
of this relation if we take the functional derivative of the energy functional 
E[v] with respect to the potential v. 



Sv(r) Sv(r) Sv(r) 5v(r) 



= (tf|n(r)|tf) = n(r) (44) 

where we used the Schrodinger equation H v \^f) = £?[u]|\E') and the nor- 
malization condition (^l 1 ?) = 1. Note that the equation above is nothing 
but a functional generalization of the well-known Hellmann-Feynman the- 
orem [7]. Now we can go to the density as our basic variable by defining a 
Legendre transform 

F[n]=E[v}- f d 3 rn(r)v(r) = (*[«]|f + W\V[v]) (45) 

where v must now be regarded as a functional of n. The uniqueness of such 
a mapping is garanteed by the Hohenberg-Kohn theorem [5]. The set of 
densities for which the functional F[n] is defined is the set of so-called v- 
representable densities. These are ground state densities for a Hamiltonian 
with external potential v. The question which constraints one has to put 
on a density to make sure that it is u-representable is known as the v- 
representability problem. We postpone a discussion of these matters to 
later sections. From SE/Sv = n it follows immediately that 

SF 

= -v(t) (46) 



5n(r) 



This is our first basic relation. In order to derive the Kohn-Sham equations 
we define the following energy functional for a system of noninteracting 
particles with external potential v s and with ground state wave function 

l*H>: 

E a [v a ] = ($[v a ]\f+V.\*[v.]) (47) 
with Legendre transform 

F s [n] = E[v s ] - J d 3 rn(r)v s (r) = ($[vJ|T|$[i>J) (48) 



and derivatives 



5^ = " W (49 > 
55) " -"• <r) (50) 



We see that F s [n] in Eq.(48) is the kinetic energy of a noninteracting system 
with potential v s and density n. For this reason the functional F s is usually 



denoted by T s . In the following we will adopt this notation. Finally we 
define the exchange-correlation functional E xc [n] by the equation 

F[n] = T a [n] + \ J d 3 rd 3 r'n(r)n(r')w(\r - r'|) + E xc [n] (51) 

This equation assumes that the functionals F[n] and T s [n] are defined on 
the same domain of densities. We thus assume that for a given ground 
state density of an interacting system there is a noninteracting system with 
the same density. In other words, we assume that the interacting density 
is noninteracting-w-representable. If we differentiate Eq. (51) with respect 
to the density n we obtain 

v s (r) = v(r) + j d 3 r'n(r')w(\r - r'|) + v xc (r) (52) 

where 

«~M = ( 53 ) 

defines the exchange-correlation potential. Now the state |$[u s ]) is a ground 
state for a system of noninteracting particles, and can therefore be written 
as an antisymmetrized product of single-particle orbitals <f%{v). If we now 
collect our results we see that we have converted the ground state problem 
into the following set of equations 

N i r r 

B H - E~2 / d 3 ^t(r)VVi(r)+ / d 3 rn(r)v(r) 
i=i J J 

+ i J d 3 rd 3 r / n(r)n(r / )w(\r-r / \) + E xc [n] (54) 
(_I V 2 + v(r) + J d 3 r'n{v')w{\v - v'\) + v xc {v))^{v) = e iVi (r) (55) 

N 

n(r) = ^|^(r)| 2 (56) 
i=i 

The above equations constitute the ground state Kohn-Sham equations [8]. 
These equations turn out to be of great practical use. If we can find a 
good approximation for the exchange-correlation energy, we can calculate 
the exchange-correlation potential v xc and solve the orbital equations self- 
consistcntly. The density we find in this way can then be used to calculate 
the ground state energy of the system. 

The exchange-correlation functional is often split up into an exchange func- 
tional E x and a correlation functional E c as 



E xc [n]=E x [n] + E c [n] 



(57) 



in which the exchange and correlation functional are defined by 



E x [n] = (*.[n\\W\* a [n\)- \ f d 3 rd 3 r'n(r)n(r')w(\r - r'\) (58) 

E c [n] = {*[n]\H v mn]) - (*M\Hv\*a[n]) (59) 

In this equation |<5> s [n]) is the Kohn-Sham wavefunction and the true 

ground state wavefunction of the interacting system with density n. Since 
the Kohn-Sham wavefunction is not a ground state wavefunction of the true 
system we sec immediately from the variational principle that E c < 0. We 
see also that the form of the exchange functional in terms of the Kohn-Sham 
orbitals is identical to that of the exchange energy within the well-known 
Hartree-Fock approximation. However, since the Kohn-Sham and Hartree- 
Fock orbitals differ, the value of E x [n] is not equal to the Hartree-Fock 
exchange. We finally remark that splitting up E xc into an exchange and 
a correlation part has several disadvantages. First of all, this splitting has 
only meaning if the ground state of the true or Kohn-Sham system is non- 
degenerate. We will later see that -E xc [n] is a well-defined functional even 
for degenerate ground states, but that exchange and correlation separately 
are ill-defined in that case. Secondly, there are many cases, notably molec- 
ular dissociation cases, where the exchange-only theory is a bad starting 
point for the treatment of correlation effects and for which it is much easier 
to find good approximations for the combined exchange-correlation func- 
tional. 

In this section we gave a derivation of the Kohn-Sham equations. However, 
in this derivation we made a number of assumptions. We assumed that 
for every ground state density of an interacting system there is a noninter- 
acting system which has the same density in its ground state. Secondly, 
we assumed that the density functional were differentiablc. This assumes 
that the values of the functionals change smoothly with changes in the 
density. In the following sections we will investigate to which extent these 
assumptions are justified. 

6. Definition of the functional derivative 

Let us start by defining what we mean with a functional derivative. The 
derivative we will talk about is what in the mathematical literature [9, 10, 
11, 12, 13, 14] is refered to as a Gateaux derivative. Let G : B — > 1Z be 
a functional from a normed function space B to the real numbers 1Z. If 
for every h e B there exist a continuous linear functional SG/Sf : B — > 1Z 
defined by: 

g[^lim G [ / + ^-^ (60) 



then SG/Sf is called the Gateaux derivative in / e B. If the limit exists 
but the resulting functional of h is not linear or continuous then this limit 
is called a Gateaux variation. Note that the definition of the Gateaux 
derivative is very similar to the definition of the directional derivative in 
vector calculus where the linear functional corresponds to the inner product 
of the gradient vector with the vector which specifies the direction of the 
differentiation. Often the linear functional SG/Sf can be written in the 
form 

jjVA = J d\g{v)h{v) (61) 
If this is the case we write 

SG 

9ir) - m (62) 

which we will call the functional derivative of G. Note further that although 
SG/Sf specifies a linear functional when acting on h, the function g(r) 
depends in general on / in a nonlinear way. We also see that the set of 
functional derivatives on a function space B is equal to the set of continuous 
linear functional on that space. This set is called the dual space of B and 
denoted B* . For instance, the dual space of B = L 1 n L 3 is known to be the 
space B* = + L°° . From this we see that if the derivative of a density 
functional defined on the set of densities L 1 n L 3 exists then its derivative 
is in the set L§ 

Let us now mention a straightforward consequence of the definition of the 
Gateaux derivative. Suppose that the functional G which we assume to 
be Gateaux differentiable, has a minimum at /o, i.e. G\f\ > G[fo] for all 
/ € B. Then the function 

g(e) = G[f + eh] > G[f ] (63) 

has a minimum in e = and therefore the derivative of g(e) in e = 
vanishes. Thus 

o -§<«)- ft *^-g(JWN <«' 

Therefore a necessary condition for a differentiable functional for having a 
minimum at /q is that its Gateaux derivative vanishes at / . We further 
prove one other fact that we will use later. If G is a convex functional, i.e. 

G[A/ + (1 - A)/i] < AG[/ ] + (1 - X)G[h] (65) 

for < A < 1 then the function 



.g(A) = G[f + A(/i - /„)] 



(66) 



is a continuous function on the interval [0, 1]. To show this we will show 
that g(X) is convex. Take Ai, A 2 and e from the interval [0, 1]. Then 



0(eAi + (1 - e)A 2 ) - G[e(f + Ai(/i - /o)) + (1 - e)(/ + A 2 (/i - /o))] 
< eff (Ai) + (1 - e) 5 (A 2 ) (67) 

Therefore 5(A) is a convex real function on the interval [0, 1] and hence 
continuous. The fact that convexity implies continuity is true for functions 
on the real axis, but this does not extend to infinite-dimensional spaces. 
Let us finally make some remarks on higher order derivatives. If g(e) = 
G[fo + eh] defines a n-fold differentiable function of e then we define the 
n-th Gateaux variation of G as 

Pft-.»i-tl«» m 

where S n G/Sf n has now n arguments h. If this expression defines a multi- 
linear continuous functional then we call this the nth-order Gateaux deriva- 
tive of G. For more details on this point we refer to reference [9]. Again, if 
this expression can be written in the form 

S n G f 

-[h,...,h] = J d 3 r 1 ...d 3 r n g(r 1 ...r n )h(r 1 )...h(r n ) (69) 

Sf(r 1 )...Sf(r n r 9iri --- rn) m 
the nth-order functional derivative of G. Sofar our discussion has been 
rather abstract. Let us therefore apply the definition and calculate some 
functional derivatives. 



Sf 

then we call 



6 n G 



7. Static linear response of the Schrodinger equation 

We now consider the effect of small changes in the external field on the 
expectation values of physical observables. This is exactly what is studied 
in most experimental situations where one switches on and off an external 
field and studies how the system reacts to this. We will here study a more 
specific case in which we look at static changes in the external potential 
and their accompanying changes in the ground state expectation values. By 
investigating this problem we will learn how to take functional derivatives 
and how the existence of these derivatives is related to the existence of the 
linear response function. 

Suppose that we have solved the following ground state problem: 



H\9 ) = (f +V + W)\^ ) = E \ * > 



(71) 



where |W ) is a nondegenerate ground state \*f?[v]) of Hamiltonian H with 
external potential V. From the wavefunction we can then calculate the 
expectation value of any operator O which is a well-defined functional of 
the external potential 

0[v] = (Mv}\OMv]) (72) 

Let us now calculate the functional derivative SO/Sv at a given potential 
v. According to our definition in the previous section we have to calculate 
the quantity 

SO.. . 0[v + eSv] - 0[v] 

— \5v] = hm — L — (73) 

6v e^o e 

To evaluate this limit we have to calculate 0[v + e<5u] which we will do using 
static perturbation theory. We therefore make a slight change 



eSV 



= e J d 3 rn(r)5v(r) (74) 



in the external potential of Hamiltonian H, i.e. we change the potential 
to V + eSV. The new ground state wavefunction which we will denote by 
|^(e)) satisfies 

(H + eSV)We)) = (75) 
We will solve this equation to first order in e with the condition |^(0)) = 
l^o)- We note that the solution of Eq.(75) is not unique, because if |^(e)) 
is a solution then also |$(e)) = e ie ^\^/(e)) is a solution, where 9(e) is an 
arbitrary function of e. If we choose 9(0) = then |4>(e)) also satisfies the 
condition |$(0)) = |W ). The arbitrariness of the phase factor obviously 
does not affect the value of any expectation value, i.e. 

0(e) = <*( c )|0|*(e)> = me)\0\m) (76) 

which is a unique function of e. However, it affects the appearance of our 
first order expansion, which for both functions looks like 

|*(e)) = |* ) + e|v]/'(0)) + O(e 2 ) (77) 
|$(e)) = |*o> + e(|*'(O)> + z0'(O)|vI/ o )) + O(e 2 ) (78) 

where |^'(0)) and 9'(0) are the first order derivatives of |^(e)) and 9(e) in 
e = 0. We see that |W(e)) and |$(e)) differ in first order by the amount 
ie9'(0)\ , i> n ), i.e. by an imaginary number times the unperturbed ground 
state. This is exactly the freedom we will find in our expansion when we 
try to obtain the first order change in the wavefunction. 
Let us expand the wavefunction and energy in Eq.(75) to first order in e. 
We then obtain 



(H - E )\*'(P)) - (E'(0) - SV)\^ ) (79) 



where E'(0) is the first order derivative of E(e) in e = 0. To solve this 
equation we expand |W'(0)) in an orthonormal set of eigenstates of the 
unperturbed Hamiltonian H as 



|*'(0)> = (80) 



?=o 



If the unperturbed Hamiltonian has a continuous spectrum then for the 
corresponding energy eigenstates the summation in this equation has to be 
replaced by an integration. If we insert the expansion (80) into Eq.(79) for 
|^'(0)) we find the equation 

oo 

- £ )|*i) - (E'(0) - SV)\* ) (81) 

i=0 

where the energies E; t for i > arc the cigcncncrgics of the excited states 
of the unperturbed Hamiltonian. If we multiply this equation from the left 
with ('Jo | w e obtain for the change in energy 

E?(Q) = {*o\SV\* Q )= J d 3 rn Q (r)Sv(r) (82) 

where no is the density of the unperturbed system. We have therefore 
shown that 



E'(0) = lim E[v + e8v]-E[v] = J 



(83) 



As has been derived before we find 5E/5v(r) = n (r). If we now multiply 
Eq.(81) from the left with (* fc | for k > we find 

Ci = _(M£M (forfc>0) (84) 

Note that these coefficients are well-defined because E k > E , since we are 
dealing with an isolated nondegenerate ground state. We therefore find for 
the first order change in the wavefunction 

i«w-«it.)-s '''><T' l ff , ' > w 

The coefficient Co remains undetermined by these equations. However, from 
the requirement that the perturbed system remains normalized we find 

= ^|^(0) = (*'(0)|*o> + (*o|*'(0)) =4 + c (86) 



and hence c must be purely imaginary. This is exactly the kind of freedom 
in the first order wave function that we noted before. However, as we have 
seen this is related to choice of phase and does not affect the expectation 
values. We now want to calculate the first order change in an arbitrary 
expectation value due to a change eSv of the external potential. This change 
is given by 

O'( ) = lim + -OM = (V(0)|6|* ) + (*o|6|*'(0)) (87) 

e— >0 e 

where O'(0) is the derivative of O(e) in e = and where we assumed that 
the operator O does not depend on V (as for instance the Hamiltonian 
does). From Eq.(85) we obtain 

O'(0) = (c5 + c )(*o|0|*o> 

^ (y \6V\y k )(* k \d\H> ) + (*o|0|* fc )<* fc |<SV1*o) 



fe=i 



E k — E 



From the normalization condition Eq. (86) we see that the first term on the 
right hand side vanishes and we can rewrite the equation as 

0>W= S £[6v] = Jd*r^- ) 5v(r) (89) 



where 



SO ^<*o|0|*fc><*fc|n(r)|*o> +&c _ 



Sv(r) ^ E k - E 

For a system with a nondegenerate ground state we therefore have obtained 
an explicit expression for the functional derivative SO/5v of the expectation 
value 0[v], evaluated at potential v(r). We will use this formula later when 
we are studying functional derivatives with respect to the density. 

8. Invertability of the density response function 

In the previous section we obtained an expression for the functional deriva- 
tive of an arbitrary expectation value. We now make a special choice for 
the operator O and we choose O = n(r'). In that case we obtain from 
Eq.(90): 

y(r ' ,fa(r')_ (tt |n(rQ|tt fc )(tt fc Kr)|tto) 

X[T > r) - Sv(r) - - Ek -E - ' - + c ' c - W 



where x(r',r) is the static density response function. From this expression 
we see immediately that \ is rea l an d symmetric, i.e. x(r,r') = x( r ', r )- 



The density response function relates first order changes in the potential to 
first order changes in the density according to 

Sn(r) = J d 3 r' X (r,r')Sv{r') (92) 

If we change the potential by a constant then the density does not change. 
We therefore must have 



dV x (r,r') = (93) 

which can also be checked immediately from Eq.(91). Because of this prop- 
erty we also see that the induced density Sn(r) automatically integrates 
to zero. The density response function has therefore indeed the physical 
properties that we expect. On the basis of the Hohcnberg-Kohn theorem 
we may further expect that the only potential variation that yields a zero 
density variation is the constant potential, i.e. if 

= J d 3 r' X (r,r')Sv(r') (94) 

then 8v(r) = C. This readily proven from the properties of the response 
function. Suppose Eq.(94) is true for some Sv(r'). Then we obtain by 
integration with 6v(r): 

/°° I |2 
d 3 rd 3 r'5v{v) x (v, t>)8v{t>) = -2 Y, (95) 
fe=i k 

where we defined a k by 

a k = J d 3 r(H> k \n(r)\H> )5v(r) (96) 

Since E k — E > equation (95) can only be true if a k — for all A; > 1. 
But this implies that 

oo oo „ 

Y a ^k)^Y / d 3 r\* k )(* k \h(r)\* Q )6v(r) 

k=l k=\ J 

= J d 3 rSv(r)(l-\^ )(^ a \)n(r)\^ a ) 

= J d 3 rSv(r) (n(r) - n (r)) |* ) (97) 

where no(r) = (^ , |«( r )| v I'o) is the ground state density. Written in first 
quantization this implies that 



= 



^A«(r,)|*o) = (98) 



i=l 



where N is the number of electrons and Av the potential 
Av(r) = 6v(r) f d 3 rrio(r)5v(r) 




(99) 



Now Eq.(98) implies that Av(r) = which together with Eq.(99) implies 
that Sv(r) — C. We have therefore shown that only constant potentials 
yield a zero-density variation, and therefore the density response function 
is invertible up to a constant. One should, however, be careful with what 
one means with the inverse response function. The response function de- 
fines a mapping \ '■ $V ~~ * SA from the set of potential variations from a 
nondegenerate ground state, which we call 5V and is a subset of L s + L°° , 
to the set of first order densities variations, which we call 5 A, that are 
produced by it. We have just shown that the inverse x 1 '■ SA — > SV is 
well-defined modulo a constant function. However, there are density vari- 
ations that can never be produced by a potential variation and which are 
therefore not in the set 5 A. An example of such a density variation is one 
which is identically zero on some finite volume. 

From our analysis we can further derive another property of the static re- 
sponse function. Since x(r,r') is a real Hcrmitian integral kernel it has an 
orthonormal set of eigenfunctions which can be chosen to be real. Let /(r) 
be such an cigenfunction, i.e. 



Multiplication from the right with /(r) and subsequent integration yields 



We therefore find that A < 0. However, we already know that A — is only 
possible if /(r) is constant. We have therefore obtained the result that if 
a non-zero density variation is proportional to the potential that generates 
it, i.e. 5n = X5v, then the constant of proportionality A is negative. This is 
exactly what one would expect on the basis of physical considerations. In 
actual calculations one indeed finds that the eigenvalues of x are negative. 
Moreover it is found that there is a infinite number of negative eigenvalues 
arbitrarily close to zero, which causes considerable numerical difficulties 
when one tries to obtain the potential variation that is responsible for a 
given density variation. We finally note the invcrtability proof for the static 
response function can be extended to the time-dependent case. For a recent 
review we refer to [15]. 




(100) 




(101) 



where 




(102) 



9. Functional derivatives and f-representability 



We are now ready to tackle the question how to calculate functional deriva- 
tives with respect to the density. By the Hohcnbcrg-Kohn theorem every 
density associated with a nondegenerate ground state uniquely determines 
that ground state and the external potential that produced it. Because 
the density n determines the external potential v and vice versa we can 
parametrize the ground state either by the density or the external poten- 
tial. The same is, of course, true for the expectation value of any operator 
O that we calculate from the ground state. We will therefore write 

0[n] = (tf [n]|6|tf[n]) = (*H|0|*H) = 0[v] (103) 

Since we know how to calculate the functional derivative with respect to 
v the functional derivative with respect to n seems rather straightforward. 
By the chain rule for differentiation we obtain 

50 ._ f jsj 50 Mr') _ f «, 50 



5n(r) 



f rfV 6 ° 6v{r) = f d\' — Y-Hv' r) (104) 
J 5v{v>) Sn(r) J 5v{v') X ( ' j [ ' 



where we used that the density response function has an inverse (modulo a 
constant). However, great care must be taken when deriving this equation. 
One needs to consider under which conditions the use of this chain rule 
is applicable. Let us therefore start with the definition of the functional 
derivative and try to calculate 

6 ^ [ Sn] = lua 0[no + e6n] -° [no] (105) 

where no is the density of a given nondegenerate ground state with external 
potential vo- Now we run into the following problem. In order to make the 
value of O[no + eSn] well-defined we have to make sure that no + eSn E 
A, i.e. we have to make sure that this perturbed density belongs to a 
nondegenerate ground state. Whether or not this is true obviously depends 
on the choice we make for Sn. A natural choice would be to take Sn E SA, 
i.e. we take a Sn from the set of first order density variations produced 
by some potential variations Sv £ SV. In this way we know that for every 
such Sn there is a potential v + eSv that generates, to first order in e, the 
required density. However, this still does not imply that n + eSn G A. 
This will only be so if we can find an additional potential that would make 
the higher order terms in e disappear and it is not clear how to prove the 
existence of such a potential. 

Nevertheless, this idea is sufficient to make a useful statement about the 
functional derivative. Let us therefore take some Sn € SA. Then, by the 



invertability of the density response function, there is a unique Sv modulo 
a constant, such that 



Sn(r) = J <Pr'xo(r,r')6v(r') (106) 

where %o is the static density response function belonging to the nondegen- 
erate ground state |^o) with density n . Consider now the set of external 
potentials vq + eSv. For each e we can solve the Schrodingcr equation and 
obtain a ground state density n e (r), which for e small enough will corre- 
spond to a nondegenerate ground state, i.e. n e (r) e A. By construction 
n e — n will, to first order in e, be equal to the chosen Sn we started with, 
i.e. n £ can be written as 

n e (r) = n (r) + e6n(r) + m e (r) e A (107) 

where the term m e (r) satisfies 

lim ^ = (108) 

Since n e is the density of a nondegenerate ground state the expectation 
value 0[n £ ] is now well-defined and we can calculate 

0[n t ] - O[no] _ O[v + eSv] - O[v ] _ f , 3 <50 



lim — ^ ^ = lim — = / d 4 r— -Sv(r) (109) 

£^o e e^o e J 5v(r) 

We have therefore found a parametrized set of densities n e £ A for which 
the limit above is well-defined. Moreover, in this equation Sv is uniquely 
defined by the Sn G 5A we started with by Sv = Xo 1 ^, so that we can 
write 

Um O[n + eSn + mt }-O[n ] = t ^SO 
e J 8n(v) 

where the functional derivative SO/5n{v) is given (modulo a constant) by 

Note that this derivative is defined up to a constant, since Sn(r) integrates 
to zero. The functional derivative SO/Sn[5n] in Eq.(llO), regarded as a 
linear functional acting on Sn, is of course independent of this constant. 
We have therefore a result that is in accordance with our naive expectation 
of Eq.(104). The functional derivative in Eq.(lll) is well-defined in terms 
of the properties of the unperturbed system and independent of the choice 
for Sn, and therefore independent of the parametrized path n £ e A that 



we used to approach n . Moreover, n f approaches the straight path n + 
e5n arbitrarily closely for e — > and this is exactly the limit that we are 
interested in. We will therefore call the SO/Sn of Eq.(lll) the functional 
derivative of 0[n\. Let us check that Eq.(lll) gives us back some known 
results. Let us take O = T + W so that 0[n] = Fhk[ti]- If we insert this 
operator in Eq.(lll) we see that we have to calculate 

(*o|T + W|* fc ) = (* |#o - Vol**) 

= E k S k0 - J d 3 rv (r)(^ \h(r)\^ k ) (112) 

where Hq = T + Vq + W . If we insert this matrix element into Eq.(lll) we 
obtain: 



SF HK [n] 
Sn(r) 



d?r'd?r"v {T")x {T" ,t')xo\t' ,t) = -v (r) (113) 



This is exactly the result that we derived earlier using Legendre transforms. 
We can, of course, also do this derivation for a noninteracting system where 
W — and therefore O = T, in which case we obtain 

,W - -«a[no](r) (114) 



5n(r) 



where v s [no] is the potential that in a noninteracting system generates 
density n a , i.e. the Kohn-Sham potential corresponding to density n . 



10. The Hohenberg-Kohn theorem for degenerate ground 
states 

Until now we considered only nondegenerate ground states. These ground 
states have the simplifying feature that they are determined uniquely (up to 
a phase factor) by the external potential and therefore any expectation value 
calculated from the ground state wavefunction is a well-defined functional 
of the external potential v. However, degenerate ground states do occur 
(for instance in open shell atoms) and in that case the external potential 
does not generate a unique ground state but a linearly independent set of 
q different ground states. The expectation value of any operator (except 
the total energy) will then depend on which ground state out of the ground 
state manifold we choose to compute the expectation value from. This is 
in particular true for the density operator, and therefore different ground 
states out of the ground state multiplct will yield different densities 
n,. If we want to consider degenerate ground states, then the density is 
no longer a unique functional of the external potential. However, we will 



show that the inverse mapping n — > v is well-defined, i.e. every ground 
state density determines uniquely the external potential that generated it. 
We will first generalize this statement somewhat. Instead of pure-state 
densities, which come from an eigenstate of the Hamiltonian H, we will 
consider ensemble densities. To define this concept we first introduce for a 
g-fold degenerate ground state {|^»), i = 1 ... 9} the density matrix 

D = 5^A i |* i )<* i | E A * = 1 (0<A,<1) (115) 

i— 1 i 

where the ground state wavefunctions are chosen to be orthonormal. 
We then define the ground state expectation value of operator O as 

(0)=TrDO (116) 

where the trace operation for an arbitrary operator is defined as 

00 

Tri = ^(0> 4 |i|<f> 4 ) (117) 

»=i 

where {|<&i)} is an arbitrary complete set of states. The trace is independent 
of the complete set that we choose, as follows simply by insertion of a 
different complete set {|*»)} 

00 00 

Tri = J2^\ A \^)= E^l 1 !^)^!^) 

1 — 1 i,j= 1 

= f;<* J .|$ < x$ i |i|* J .) = f;<* j |i|* j ) (us) 

If we choose the complete set to be the set of eigenstates of Hamiltonian H 
then we find 

00 q 

TvbO = J2(*i\DO\*i) = ^i\d\*i) (119) 



This defines the expectation value of an observable O in an ensemble de- 
scribed by density matrix D. For the particular case of the density operator 
we have 



q 1 
n(r) = TrDn(r) = ]T i\h(r)\H> J = £ A,n 4 (r) (120) 

i=l i=l 



We will denote densities n(r) of this type, which are obtained from an 
orthonormal set of ground states {l^j)^ = 1 . . . q} corresponding to a po- 
tential v, as ensemble w-representable densities, or for short E-V-densities. 
We further denote the set of all E-V-densities generated by potentials in 
Li + L°°as B. A density will be called a pure state ii-representable den- 
sity or for short PS-V-density if it can be written as n(r) = (W|n(r)|W) , 
where |W) is a ground state. Obviously PS-V-densities are special cases of 
E-V-densities and the set of PS-V-densities is therefore a subset of the set 
of E-V-densities. 

We are now ready to formulate the Hohenberg-Kohn theorem for degen- 
erate ground states: Every E-V-density determines the external potential 
that generated it up to an arbitrary constant. Let us make this statement a 
bit more specific. Suppose that D\ and D 2 are density matrices belonging 
to ground state ensembles for potentials v\ and v 2 respectively, with corre- 
sponding ensemble densities m and n 2 . If v\ ^ v 2 + C with C a constant, 
then m ^ n 2 . The proof is analogous to the proof of the nondegenerate 
case. 

Suppose vi generates the ground state multiplct A± = i = 1 ■ ■ • qi} 

and v 2 generates the ground state multiplct A 2 = {\^fi),i = 1 . . . q 2 }. All 
the wavefunctions within these multiplets may, without loss of generality, 
be chosen orthonormal. Then none of the wavefunctions in the sets A\ and 
A 2 are equal. This follows from the same argument as used in the proof of 
the Hohenberg-Kohn theorem for the nondegenerate case. In particular, as 
the sets A\ and A 2 are only defined to within a unitary transformation no 
in A 2 is a linear combination of the in A\. This then implies that 
two ground state ensemble density matrices constructed from the ground 
states in A\ and A 2 are different 

91 92 

L\ = A<|*i><*i| + X>|*i)<*i| - D 2 (121) 

i=l i=l 

where ^2i \ = ^2 Hi = 1. This follows, for instance, by taking the inner 
product on both sides with |^ m ) as the \*S>i) are not linear combinations 
of the We have thus established that the sets of ground state density 
matrices for the two different potentials v\ and v 2 are disjoint. We now have 
to prove that the density matrices in these sets yield different densities. 
If Hi = f + Vi + W and H 2 = f + V 2 + W then 

TvDiH 2 > TrD 2 H 2 (122) 

This follows directly from 



91 91 

i=i »=i 



= J2\E[v 2 ] = E[v 2 ] = J2^i\H 2 \^i) = Tr D 2 H 2 (123) 

i=l i=l 

We can now show that D\ and D 2 yield different densities. We proceed 
again by reductio ad absurdum. We have 

E[ Vl ] = TrD 1 H 1 = TiDi(H 2 + Vi-V 2 ) 

= Tr DiH 2 + J m(r)(v\(r) — v 2 (r))dr 

> TrD 2 H 2 + [ m(r)(«i(r) - v 2 (r))dr 



Likewise we have 



= £7[v2] + y"ni(r)(vi(r)-«2(r))dr (124) 



E[v 2 ] > E[ Vl ] + J n 2 (r)(v 2 (r) - Vl (r))dr (125) 

which added to the last inequality leads to 

J d 3 r(n 2 (r) - m(r))(« 2 (r) - Wl (r)) < (126) 

This leads again to the contradiction < if we assume that n\ = n 2 . 
Therefore D\ and D 2 must give different densities, which proves the theo- 
rem. 

Within the set of ensemble ground state density matrices corresponding 
to the same potential however, two different density matrices can yield the 
same density. The simplest example is the hydrogen atom with the degener- 
ate ground state wavefunctions for an electron with spin up and spin down, 
i.e. W| = (Is) a and ^ = (ls)0 where a and (3 are the spin wavefunc- 
tions. Both wavefunctions obviously have the same total density. Another 
example is provided by the the two degenerate ground state wavefunctions 
(\s) 2 2p + and (ls) 2 2p~ of the noninteracting lithium atom, where 2p ± are 
p-orbitals with angular momentum quantum numbers I = ±1. Both wave- 
functions lead to the same density 

n(r) = 2|^ ls (r)| 2 + |^ 2p+ (r)| 2 = 2|^ ls (r)| 2 + |^ 2p -(r)| 2 (127) 

We have therefore constructed two degenerate ground states with the same 
density. As a consequence the ground state expectation value of a given 
operator may no longer be considered as a functional of the density (take, 
for instance, the expectation values of the spin and angular momentum in 
our examples). However, if two different ground state density matrices have 



the same density then also the energy TrDH for those different density ma- 
trices is the same. For every E-V-density n we can therefore unambiguously 
define the ensemble version of the Fhk functional as [16] 

FehkH = TrD[n](f + W) (128) 

where D[n] is any of the ground state ensemble density matrices corre- 
sponding to n. We can now define an extension of the energy functional 
E v to the set of E-V-densities 

E v [n] = J n(r)v(r)dr + F EHK [n] = Tr D[n]H (129) 

Similarly as for Fhk we easily can prove 

E[v] = inf |y n(r)v(r)dr + F EHK [n]j (130) 

The functional Fehk is an extension of Fhk since 

FehkH = F HK [n] if neA (131) 

This follows directly from the fact that for a non-degenerate ground state 
\^f[n]) corresponding to n we have D[n] = \^f[n])(^f[n]\, so 

F EHK [n]=TrD[n](f+W) = (*[n]|f + W|*[n]) = F HK [n] (132) 

We can furthermore prove that Fehk is convex by the same proof as for 
Fhk- Nothing is however known about the convexity of the set of E-V- 
densities B itself which constitute the domain of Fehk- Let us collect our 
results in the form of a theorem 

Theorem 2 (Hohenberg-Kohn) The E-V-density n specifies the exter- 
nal potential up to a constant. Moreover, 

1. The ground state energy of a system with external potential v can be 
obtained from 

E[v] = inf I / n(r)v(r)dr + FehkM 

neB [J 

where F EHK [n] = TrD[n](T + W) and D[n] is a ground state density 
matrix with TrD[n]h(r) = n(r). 

2- Fehk is convex 



As we will now demonstrate the subset of PS- V-densities of B is not convex. 
More precisely, we will now show that there are E- V-densities which are 
not PS- V-densities [1, 17]. As any E-V-density is a convex combination 
of PS-V densities this then demonstrates the non-convexity of the set of 
PS- V-densities. 

Consider an atom with total angular momentum quantum number L > 
which has a 2L + 1-degenerate ground state. The external potential v is 
the spherically symmetric Coulomb potential of the atomic nucleus and 
the degeneracy is a result of rotational invariancc of the Hamiltonian of 
the system. The ground state wavefunctions then transform among one 
another according to a 2L + 1-dimcnsional unitary representation of the 
rotation group. We assume that there is no accidental degeneracy If we 
denote the ground state wave functions by {^[rij]} = = 1 . . . 2L + 1} 

and the corresponding electron densities by rij then the following convex 
combination 



is invariant under all rotations and therefore spherically symmetric. How- 
ever the rii are not spherical. In fact, there is no linear combination of the 
ground states that leads to a spherically symmetric density. As the nj 
is obtained from \^j), which by a unitary transformation can be obtained 
from any other and the external potential is invariant under rotations 
we find that 



for all < i,j < 2L + 1. Let us now suppose that n is a ground state 
density obtained from a pure state wavefunction |^[n]). This wavefunction 
is not a linear combination of the \*S?i) otherwise n would not be spherically 
symmetric. We then find 




2L+1 



(133) 




(134) 





2L+1 



(135) 



This then gives 




2L+1 



(136) 



But we also know that F EHK is convex on the set of E-V-densities. This 
leads to a contradiction and hence we must conclude that n can not be a 
pure state density of any potential. The density n is, however, a convex 
combination of ground state densities corresponding to the same external 
potential and therefore, by definition, an ensemble u-representable density. 
We therefore have constructed an E-V-dcnsity which is not a PS-V-density. 
For an explicit numerical example of such a density we refer to the work of 
Aryasetiawan and Stott [18]. 

11. Linear response for degenerate ground states and 
the differentiability of F E hk 

In the previous sections we learned to take functional derivatives on the 
basis of linear response theory. We concentrated on nondegenerate ground 
states for which the response functions are well-defined. We will now use 
response theory for degenerate ground states to study the differentiability 
of the functional Fehk- 

It is clear that in the case of degeneracy the expectation value of all ob- 
servables, except the energy, depends on which particular ground state we 
choose to calculate the expectation value from. This poses a clear difficulty 
for the definition of general density functionals and their functional deriva- 
tives. One may, however, argue that if a potential vq leads to a degenerate 
ground state then we can always find an arbitrarily small perturbation of 
the potential eSv that lifts the degeneracy and therefore the expectation 
value 

O[v + eSv] = (^[v + eSv}\6\^[v + eSv}} (137) 

of an observable described by an operator O does exist for any e > 0, where 
^[uo + eSv]) is the ground state wavefunction of the perturbed system. 
Therefore also the limit 

O [v , Sv] = lim O [v + eSv] (138) 

e— >0 

is defined, although it will in general depend on the potential variation Sv. 
One can therefore define the functional derivative as follows 

O>[5v] = Km O[V0 + e6v] -° [V0 ' Sv] (139) 

However, we will see that these derivatives depend in a nonlinear way on Sv 
and are therefore only defined as Gateaux variations rather than Gateaux 
derivatives. We shall especially be interested in the case where O is the 
density operator. For the energy the limit in Eq.(138) is independent of Sv 
and we may wonder if a Gateaux derivative still exists. Let us investigate 



this in more detail. Suppose we have an Hamiltonian H with external 
potential vo which has a g-fold degenerate ground state. Let us now apply 
a perturbation eSv which lifts the degeneracy. Then there arc q cigcnstatcs 
{|*fc( e ))>fe = 1 • • • l} with energies E k (e) of the perturbed Hamiltonian 
that are continuously connected to degenerate ground states |Wfe(0)) of the 
unperturbed Hamiltonian H , i.e. 

lim |* fc (e)> = |* fc (0)) (140) 

£—►0 

with 

E k (0) = (*fc(0)|# o |*fe(0)) = E (141) 

Note that since the states |\l/fc(e)) are orthonormal, also the limiting states 
states |^fe(0)) in the degenerate ground state multiplet are orthonormal. 
Which particular ground states of the ground state manifold are reached 
by the e — > limit depends obviously on the form of the perturbation 5v. 
In which way it depends on Sv we will investigate now. The ground state 
wavefunctions |^fc(e)) satisfy the equation 

(H + e8V+j6U)\* k (e)) = E k (e)\* k (e)) (142) 

where, for the later discussion, we also added a one-body potential SU 
which is of second order in e. It turns out that this term influences the 
first order density response. We note, like in the nondegenerate case, that 
if |Wfe(e)) is a solution to this equation then also \Q k (e)) — e l6k ^\^> k (e)) is 
a solution, where 9 k (e) is an arbitrary function of e. This freedom will, as 
in the nondegenerate case, not affect any of the expectation values. If we 
expand the Schrodingcr equation (142) to first order in e we obtain 

(H E )\%(0)) = (E' k (0) 8V)\9 k (0)) (143) 

where E' k (0) and |W^.(0)) are the first order derivatives of E k (e) and \^ k (e)) 
with respect to e in e = 0. In order to solve this equation for 1^^,(0)) we 
expand this quantity in the q e-connected ground states and all the other 
cigcnstatcs of H , i.e. 

oo 

\%m = j2 c i\*i) ( 144 ) 

1=1 

where = |^i(0)) for i = 1 ...q and the states > q} are eigen- 

states of H with eigencnergics Ei > E a . If we insert this expansion into 
the first order Schrodinger equation (143) we obtain the equation 

oo 

]T c^Ei - £ )|^> = (^(0) - <yv-)|* fc (0)) (145) 

i=l 



If we multiply this equation from the left with (*f>j | where 1 < j < q we 
obtain 

= E' k (0)6 jk -(9 j (0)\5V\9 k (0)) (146) 

This is the equation tells us exactly how Sv picks out the ground states in 
the degenerate multiplct. They are exactly the ones that diagonalize SV 
within the g-dimensional space of degenerate ground state functions. This 
equation also tells us that 

E' k (0) = {* k {0)\SV\* k {0)) = J d\n k {v)Sv{v) (147) 

If the state perturbed state with the lowest energy has label k = 1 then we 
see that the functional derivative of the ground state energy is given as 

SE 

= m(r) (148) 



5v(r) 



However, since n\ implicitcly depends on the perturbation Sv which picked 
out a particular ground state, this is not a proper Gateaux derivative but 
only a Gateaux variation. So even the energy functional has no functional 
derivative in the proper sense. Let us investigate the consequences for 
different expectation values. We therefore first have to determine the first 
order wavefunction. If we multiply Eq.(145) from the right with (vf^l where 
i > q we obtain the coefficient for i > q and the following expression for 

ino>»: 

l«i(0)> =£■*!».)-£ (149) 

i=\ »=9+l * ° 

However, the coefficients for 1 < i < q are not determined by Eq.(143). 
We know from the orthonormality of the states |\ffc(e)) that 

= ^^(0) = (*'i(0)|**(0)) + (*i(0)|*' fc (0)) - 4* + 4 (150) 

i.e. the coefficients form an anti-Hermitian matrix. In particular we find 
that the coefficients c\ are purely imaginary. The value of these diagonal 
elements is undetermined by the Schrodinger equation as they are related 
to the arbitrary phase ^(e) which we can choose for each |^ fc (e)). As in the 
nondegenerate case the value of any expectation value docs not depend on 
the coefficients c\ . This leaves us with the determination of the off-diagonal 
terms for i ^ k. These coefficients can be found by expanding Eq.(142) 
to second order in e 2 . This yields 



(H £ )K'(0)) - 2(^(0) - 5V)\%(0)) + (££(0) - ^|* fc (0)) (151) 



where ^'(O)) and E'£(Q) are the second order derivatives of |\I/fc(e)) and 
Ek(e) in e = 0. If we multiply this equation from the left with (W;(0)| where 
1 < I < q and use Eq.(146) and Eq.(149) then we obtain 



= 2(^(0) -^(0))cf + 4'(0)J, fc -<*,(0)|^|*fc(0)) 

Ei — Eq 



+ 2 y (152) 



i=q+l 

If we in this equation take k = I we obtain an expression for the energy to 
second order in e: 

Bfc(0) = 2 /" d 3 rd 3 r'^(r)xfc(r,r')^(r')+ / d 3 rn fe (r)<5u(r) (153) 

where is the density corresponding to |^fc(0)) and \k the following 
function 

x r> y ^(0)|n(,)|^)^|n(,Q|^(0)) +c& 

t=9+l 1 

If, on the other hand, we take I ^ k we obtain our desired equation for the 
coefficients cf: 

fc 1 (^(O)l^l^fc(O)) 

1 2 E' k (0) - El(0) 

££(0)-£?(0) ^-£ 1 J 

We see that indeed cf* = —c l k as expected. We have therefore completely 
determined the first order change in the wavefunction. Let us see what this 
implies for the observables. The first order change in the expectation value 
of an observable O is given by 

O'(0) = <* fc (0)|6|* , fc (0)) + c.c. (156) 

If we insert our expression for the first order change in the wavefunction 
we obtain 

O'(0) = J d 3 rCk(r)Sv(r) + l - J d 3 r Vk (r)Su(r) 

+ J d 3 rd 3 r'^ k {r,r')6v{r)6v{r') (157) 



where we defined the following functions 

Ei — Eq 



= £ M^»M +C ., (158) 



i=q+l 

q 



( , = v (M°M!iM + cc (159) 

&( ' } = (^oFWE^Eo) (160) 

If we choose O — n(r") then we obtain an expression for the first order 
change in the density 

Sn k (r") = Jd\ X k{v" 1 v)5v{v)+ l -Jd\ m {v" 1 v)8u{v) 

+ j rf 3 rd 3 r'^(r",r,r')<5w(r)^(r') (161) 

where the functions \k, Vk an d S,k are obtained by inserting n(r") for O in 
Eqns.(158), (159) and (160). If we denote the perturbed state with label 
k = 1 as the one with the lowest energy then 5ri\ represents the change 
in ground state density. We see that this density change does not depend 
linearly on the first order change Sv in the potential and even the second 
order change Su of the potential contributes to its value. Let us now take 
Su = and write 

<5ni(r") = J d 3 r X i(r",r)8v(r) 

d 3 rd 3 r'Zi(r", r, r')Sv(r)6v(r') (162) 



We can ask ourselves the question whether a given first order variation 8n\ 
uniquely determines the first order density change Sv. One can show from 
the Hohcnberg-Kohn theorem for degenerate states that this is indeed the 
case. If we in Eq.(126) take vi = v\ + eSv(r) where 5v is not a constant 
function we obtain 

e 2 J d 3 r5 ni (r)Sv(r) + 0(e 3 ) < (163) 
Dividing by e 2 and taking the limit e — > then yields 

d 3 rSn 1 (r)Sv(r) < (164) 



We see that the first order density change Srii can not be zero if Sv is not 
a constant. In contrast to the nondcgcneratc case Eq.(164) does not imply 
that the function \i is invertable, only the relation in Eq.(162) is invertable 
where the inverse only exists on the set of u-representable density variations 
within the set B. 

Let us now see what the linear response theory can tell us about the differ- 
entiability of F EHK . We consider again a system with a degenerate ground 
state and external potential vq. Let us take out of the set of ground state 
ensemble densities for this system a particular density tlq. We may then 
wonder whether the following limit 

SF EHK F EH K[no + eSn] - F EHK [n ] 

— \Sn = lim 165 

exists. For this limit to be well-defined we have to make sure that n + 
eSn g B and we run into the same problem as for the nondegcnerate case. 
However, in the nondegenerate case we could solve this difficulty using 
response theory. We will do the same for the degenerate case, but we will 
see that some difficulties remain. Suppose we look at the perturbed system 
with potential v + eSv in which the degeneracy is lifted. For this system 
there is a well-defined ground state density n e = n[vo + eSv]. From response 
theory we known that we can write n e as 

n e (r) = m(r) + eSn^r) + m e (r) (166) 

where 

ni(r) = lim n[v a + eSv] (167) 

= lim^^ (168) 

and where 8ni is explicitly given in Eq.(162). Now by construction the 
density n e is in B for all values of e. We can now consider the limit 

SF EHK F E Hx[ n i + tSni + to £ ] — F E hk{tii] 

-[oni\ = lim 



Sn e^O 



e 



lim f ehk[vq + eSv\ - F EHK [v , Sv] 

where F E hk[vq, Sv] is defined as in Eq.(138). The latter limit is readily 
calculated by inserting T + W = Ho — Vq for the operator O in Eqns.(158) 
and (160). We obtain for the first order change of F E hk 



SF ehk = - j d z rd z r"v {v") X i{v'\v)5v{v) 



- J d 3 rd 3 r'd 3 r"v {r")£ 1 {r",r,r')6v{r)6v{r') 

= -J d 3 rv (r)S ni (r) (170) 

and we thus obtain 

,. F EHK [m + edn^ + m^- F EHK [m] f 3 

lim — - = - d a rv (r)Sn 1 (r) (171) 

We see that this limit is linear in Sn\ and that furthermore vo is independent 
of 8n\. Therefore we can write for the functional derivative of Fehk in n\\ 

w (ni) = ^ o(r) (172) 

up to an arbitrary constant. However, we see that this derivative only exists 
for some special set of ground state densities corresponding to potential vq. 
The derivative exists for those ground state densities that correspond to 
pure states that can be obtained in the e — > limit for a perturbed system 
with potential vq + eSv. It is not clear how to take the functional derivative 
at an arbitrary ensemble density n for potential v . We saw that there 
exist E-V-densities that are not PS-V-densities. If we consider an ensemble 
corresponding to such a density and change the external potential to v +eSv 
where Sv lifts the degeneracy then for e > the ensemble will change into a 
pure state and the density will change abruptly. We must therefore conclude 
that for general E-V-densities the functional derivative of Fehk does not 
exist. This poses not only a theoretical problem but also a practical one. As 
we will discuss later, it is known from numerical investigations that there 
are ground state densities of interacting systems that are not pure state 
densities for a noninteracting system and therefore there is a clear need for 
establishing a Kohn-Sham scheme for arbitrary E-V-densities. Fortunately 
it turns out that one can define an extension of the functional Fehk to a 
larger domain of densities which can be shown to be differentiable at the 
set of all E-V-densities. This functional is the Lieb functional F^ and will 
be studied in the next section. 



12. The Levy and Lieb functionals F LL [n] and F L [n] 

The functionals Fhk and Fehk have the unfortunate mathematical dif- 
ficulty that their domains of definition A and £>, although they are well- 
defined, are difficult to characterize, i.e. it is difficult to know if a given 
density n belongs to A or B. It is therefore desirable to extend the domains 
of definition of F HK and F EHK to an easily characterizable (preferably 



convex) set of densities. This can be achieved using the constrained search 
procedure introduced by Levy [19]. We define the Levy-Lieb functional Fll 
as 

F LL [n] = inf (*|T + W\V) (173) 

where the infimum is searched over all normalized anti-symmetric iV-particle 
wave functions in iJ 1 (7?. 3Ar ) yielding density n. As shown earlier such a den- 
sity is always in the convex set S which is again a subspace of L 1 n L 3 . One 
can furthermore show, as has been done by Simon [1], that the infimum is 
always a minimum, i.e. there is always a minimizing wave function. As this 
is the first important result that we obtain for F LL we put it in the form 
of a theorem 

Theorem 3 For any ne5 there is a |*[n]) £ H 1 (1l 3N ) such that 
FllM = (*[n]|f+wn*[n]) 

Let us discuss some properties of -Fll- The functional Fll is an extension 
of the Hohenberg-Kohn functional Fhk , which is defined on A, to the larger 
set S, i.e 

F LL [n) = F HK [n] if n e A (174) 

This is readily derived. Suppose n is some ground state density correspond- 
ing to some external potential v and ground state |^[n]) then 



/ 



n(v)v(v)dv + F HK [n] = (V[n]\H\V[n]) = inf 



J n(r)v(r)dr + inf (*|f + W\V) = J n(r)v(r)dr + F LL [n] (175) 
We define a corresponding energy functional 

E v [n] = j ' n{v)v{v)dv + F LL [n] (176) 

If n is the ground state density for potential v with corresponding ground 
state wave function ^[n ] then 

E v [n] = inf > mn ]\H\Mn ]) - E v [n ] (177) 

Minimizing E v over the set 5 therefore yields the ground state density n 
corresponding to external potential v. The functional Fll has however the 
inconvenient property that it is not convex. 

Theorem 4 The functional F LL is not convex 



To show this we take the example of a previous section where we presented 
a density n which did not correspond to a ground state wavefunction. It 
was a convex combination of 2L + 1 degenerate ground state densities rii 
with corresponding ground states ^[n,]) for an external potential v. Then 
we find 



/ 



n(v)v{v)dv + F LL [n] = inf (*|JT|tf) 



1 2L+1 1 2L+1 „ 

> 2iTT E <*MI#I*N> ^ ^TTT E + / »(r)«(r)* 

i=i »=i J 

(178) 

and we find 

1 N 

F ^>^lH F M ( 179 ) 

i=l 

which proves the non-convexity of Fll- This is somewhat unfortunate as 
convexity is an important property which can be used to derive differen- 
tiability of functional. We will therefore now define a different but related 
convex functional with the same domain S. This is the Lieb functional 
defined as 

F L [n]= inf TrD(f + W) (180) 
where the infimum is searched over all TV-particle density matrices 

£> = ^A i |* i )<* i | £> = 1 Weff 1 ^) (181) 

i=l i=l 

which yield the given density n(r) = Trl)n(r) where is an orthonor- 

mal set. Also for this case one can prove the infimum to be a minimum, 
i.e. there is a minimizing density matrix. We put the result again in the 
form of a theorem 

Theorem 5 For every ne5 there is a density matrix D[n] such that 

F L [n] = TrD[n]{f + W) 

For the proof of this theorem we again refer to Lieb [1]. The functional Fl 
is an extension of Fehk to the larger set S, that is 

F L [n] = F EHK [n] if n e B (182) 



This follows directly from the fact that if n e B then there is a potential v 
which generates a ground state ensemble density matrix D[n] which yields 



n. So 



J n(r)v(r)dr + F EHK [n] = J n(r)v(r)dr + Tr D[n](f + W) 



= Tr D[n]H = inf TtDH 



J n(r)v(r)dr 



iv + F L [n] (183) 

We can again define an energy functional 

E v [n] = J n(r)v(r)dr + F L [n] (184) 

which by a similar proof as for Fll assumes its minimum at the ground 
state density corresponding to potential v. We further have the following 
relations 

F L [n] = F LL [n] if n e A (185) 

and 

F L [n] < F LL [n] if n e B and n A (186) 

The first relation follows from the fact that if the density n is a pure state 
f-rcpresentable density then the minimizing density matrix for Fl is a pure 
state density matrix. The second relation also easily follows. We take n to 
be an E-V-density which is not a PS-V density. There is a ground state 
ensemble density matrix D[n] for which we have 

J n(r)v(r)dr+F L [n] = J n(r)w(r)dr+Tr£)[n](f +W) = (^i\H\^i) (187) 

where is any of the ground states in the degenerate ground state multi- 
plet. Any wave function yielding density n can not be a linear combination 
of these ground state wave-functions otherwise n would be pure state v- 
representable. Therefore its expectation value with the Hamiltonian must 
be larger, i.e 

< inf (*|iT|*> = I n{v)v{v)dv + F LL [n] (188) 
which proves our statement. 

We will now demonstrate another important property of Fl , which is its 
convexity. 



Theorem 6 The functional F L is convex 



This is easily shown. If n = A 1 n 1 + A 2 n 2 with Ai + A 2 = 1 and < Ai, A 2 < 1 
then we have 



AiF L [m] + X 2 F L [n 2 ] = Ai _inf TrD^T+W) 
+ A 2 inf Tri5 2 (T + W) 

D2 — >ri2 

inf Tr(X 1 D 1 + X 2 D 2 )(f + W) 

Di ,£>2— >m ,n 2 

> inf TrD(f + W) = F L [n] (189) 

We therefore now have established that Fl is a convex functional on a 
convex space. This is important information which enables us to derive 
the Gateaux differentiability of the functional Fl at the set B of ensemble 
w-representable densities. We will discuss this feature of Fl in the next 
section. 

Having obtained some desirable convexity properties of Fl we try to obtain 
some analytic properties of this functional. An obvious question to ask is 
whether this functional is continuous. To be more precise, suppose that 
a series of densities n k approaches a given density n in some sense, for 
instance ||n — n k \\i — > and \\n — n fc || 3 — > for k — > 00 in L 1 n L 3 . Does 
this imply that |-Fi,[nfc] — — » 0. ? It turns out that this question is 
not easily answered. However, one can prove a weaker statement. Suppose 
n is an E-V-density corresponding to potential v of Hamiltonian H v . Then 

F L [n} + J d 3 rn(r)v(r) = TrD[n]H v < TrD[n k ]H v 

= F L [n k ] + J d?rn k {v)v{v) (190) 

and therefore 

F L [n] < F L [n k ] + J d 3 rv(r)(n k (r) - n(r)) (191) 

Further, since n k — > n in the norms on L 1 fl L 3 for each e > there is an 
integer M such that \\n — n k \\\ < e and ||n — n k \\^ < e for k > M. Now we 
can split up v € Li + L°° asn = « + w where it € Ls and w G i°° and we 
have 



/ 



d 3 rv(r)(n k (r) - n(r))\ < \\u\\ 3 ||n - n fc || 3 + ||w||oo II" - n k \\ 1 (192) 



So if ||n - rife] ||i < e/C and ||n - n fc ]|| 3 < e/C, where C = ||u||| + ||w||oo 
then 

J^M^FifaJ+e (193) 



for k sufficiently large. We therefore see that if we take e — > from above 
then -Ft [rife] approaches Fz,[ri] from above and rife — > n in the norms on 
L 1 C\L 3 . This, of course, does not imply that Fz,[n] is continuous in n. For 
that we would have to prove that |-Fi,[n] — Fl [rife] | < e rather than Fj,[n] — 
Fz,[rifc] < e - What we have proven is a weaker form of continuity, known 
as semicontinuity. Because the limit point is a lower bound, the functional 
with the property in Eq.(193) is called lower semicontinuous. This can 
also be characterized differently. If we define inf F[n m ] by infFi[n m ] = 
inf{.Fk[nfe]|fc > m} then lower semicontinuity implies 

F L [n] < lim infF L [n m ] (194) 

m^oo 

Since we required that n is an E-V-density we have proven that Fl is lower 
semicontinuous on the set of E-V-densitics. It turns out that one can prove 
that Fl is lower semicontinuous on all densities in S (see theorem 4.4 of 
Lieb [1]). However, the proof of this is not simple and we will therefore not 
try to reproduce it here. Since the result is important we present it here in 
the form of a theorem 

Theorem 7 Suppose rife and n e S and — ► n for k — > oo in the norms 
on L 1 n L 3 . Then 

F L [n] < lim infF L K] 

fe — >oo 

In other words Fl is lower semicontinuous on the set S. 

One can prove an even stronger theorem in which we only need weak con- 
vergence of the series rife. We will, however, not need that property in 
the remainder of this review. The notion of lower semicontinuity is an 
important property for convex functionals which will allow us to make sev- 
eral other useful statements about other properties of Fl- One particular 
consequence of lower semicontinuity that we will use, is that for a lower 
semicontinuous convex functional G : B — > 1Z the following set of points 

cpi(G) = {(n, r) e B x K\G[n] < r} (195) 

is closed and convex [11]. This set is called the epigraph of G and is simply 
the set of points that lie above the graph of G. The closedness property 
means that if a series of points (rife, rj.) in the epigraph of G converges to a 
point (n, r) then this point also lies in the epigraph. This also implies that 
the set of interior points of cpi(G), i.e. the set of points that lie strictly 
above G, is an open convex set (this means that for every point in this set 
there is a neighborhood that contains it). We will need this property in the 
next section. 



13. Differentiability of F L 



In this section we will prove that the Lieb functional is differentiable on the 
set of E-V-densities and nowhere else. The functional derivative at a given 
E-V-density is equal to — v where v is the external potential that generates 
the E-V-density at which we take the derivative. To prove existence of the 
derivative we use the geometric idea that if a derivative of a functional G[n] 
in a point no exists, then there is a unique tangent line that touches the 
graph of G in a point (n ,G[n ]). To discuss this in more detail we have 
to define what we mean with a tangent. The discussion is simplified by 
the fact that we are dealing with convex functionals. If G : B — > 1Z is a 
differentiable and convex functional from a normed linear space B to the 
real numbers then from the convexity property it follows that for no, n\ e B 
and < A < 1 that 



G[n + A(m - n )] - G[n ] < A(G[m] - G[n ]) (196) 
From the fact that G is differentiable we then find 

n\ i n\ l v G[n + A(ni - n )] - G[n ] SG 

G m - G n > hm = — m - n ] (197) 

A^o A on 

Therefore 



G[m] > G[n ] + -=-[m- n ] (198) 



Now 8G/dn is a continuous linear functional which will be identified with 
a tangent. This is the basis of the following definition. If for a convex 
functional G there is a continuous linear functional L : B — > 1Z such that 



G[m] > G[no] + L[m - n ] (199) 

then L is called a subgradicnt or tangent functional in no- The geometric 
idea behind this definition is illustrated in Figure 1 . 




Figure 1: The convex functional F has a unique tangent 
functional L in the point n 



For the Lieb functional Fl we will now prove the following statement: 

Theorem 8 The functional Fl has a unique tangent functional for every 
E-V -density and nowhere else. Moreover the tangent functional at an E-V 
density n can be identified with —v where v is the potential that generates 
this density. 

In order to show this we first define the energy functional 

E[v] = miTrDH v (200) 

D 

where the infimum is searched over all density matrices of the form given 
in Eq.(181). This is an extension of the previously defined functional in 
Eq.(130) to all potentials in Li + L°° . One can further show that if an 
infimum exists the the minimizing density matrix satisfies the Schrodingcr 
equation [1]. Let us first suppose that n is an E-V-density corresponding 



to potential v of Hamiltonian H v . Then 

F L [n }+ J d 3 rn (r)v(r) = TrD[n ]H v < Trt)[n}H v 

= F L [n]+ J d 3 rn{v)v{v) (201) 

and therefore 

F L [n) > F L [n Q ] - J d 3 rv(r)(n(r) - n (r)) (202) 

Therefore the functional L 

L[n] = - J d 3 rv(r)n(r) (203) 

defines a linear and continuous functional (continuity follows from Eq.(192) 
) and is therefore a tangent functional. Now we have to show its uniqueness. 
Suppose L' is a different tangent functional. Since F is defined on B = 
L 1 n L 3 the tangent functional must represent an element in the dual space 
B* = L 'z + L°° and can be written as 

L'[n] = - J d 3 r5(r)n(r) (204) 

with v e L'i + L°° and v ^ v + C. Now since we assumed that L' was also 
tangent functional we have 

F L [n] > F L [n Q ] - J d 3 rv(r)(n(r) - n (r)) (205) 

and hence 

F L [n] + J d 3 rn(r)v(r) > F L [n ] + J d 3 rn (r)v(r) (206) 

We take the infimum over all densities in S on the left hand side. Since 
we know that there is ground state density matrix corresponding to no we 
obtain from Eq.(206): 

E[v] = inf {F L [n] + [ d 3 rn(r)i(r)} > TrD[n ]Hv (207) 
nes J 

Now there are two cases, either the Hamiltonian with potential v is able to 
support a bound ground state or it is not. If it does support a ground state 
then 

Tr£»[n ]F 5 > E[v] (208) 



since D[n a ] is not a ground state density matrix for v. Together with 
Eq.(207) this immediately leads to the contradiction E[v) > E[v\. If the 
Hamiltonian with potential v does not support a bound ground state then 
there is no normalized density matrix for which the infimum of E[v] as in 
Eq.(200) is attained. We then again obtain Eq.(208) and the same contra- 
diction E[v] > E[v\. We therefore conclude that the E-V-density n has a 
unique tangent functional equal to — v where v is the potential that yields 
ground state density hq. 

Let us now suppose that no is not an E-V-dcnsity. Let us further suppose 
that there is a tangent functional at no, i.e. that for some v Eq. (205) 
and therefore also Eq.(206) are satisfied. The constrained search for Fl[uo\ 
always has a minimizing density matrix, which we call _D[n ]. With this 
density matrix then also Eq.(207) is true. Since no is not an E-V-density 
-D[no] can not be a ground state density matrix for the Hamiltonian with 
potential v. With the same arguments as before we find that Eq.(208) must 
be true which again leads to the contradiction E[v] > E[v]. Therefore there 
is no tangent functional at n if n is not an E-V-density. We have there- 
fore proven our statement, there is a unique tangent functional at every 
E-V-density and nowhere else. 

As a next step we will show that this implies that the functional Fl Gateaux 
diffcrcntiable at every E-V-density and nowhere else. For the application 
of the next theorem it is desirable to extend the domain of Fl to all of 
L 1 n L 3 . We follow Lieb [1] and define 

F T In] { T±D(f +W) if n e S , , 

The reader may seem surprised with the appearance of +00 in this def- 
inition. However, infinite values are well-defined in the theory of convex 
functionals [11] and they are usually introduced to deal in a simple way 
with domain questions. With this definition the functional Fl is a convex 
lower semicontinuous functional of the whole space L 1 n I? . We are now 
ready to introduce the following key theorem which we will use to prove 
differentiability of F L at the set of E-V-dcnsities: 

Theorem 9 Suppose F : B — > 1Z is a lower semicontinous convex func- 
tional which is finite at a convex subset V C B of a normed linear space 
B. If F has a unique continuous tangent functional L : B — > 1Z at n € V 
then F is Gateaux differentiable at n and L = SF/Sn. 

The fact that uniqueness of a tangent functional implies differentiability is 
clear from a geometric point of view. In Figure 2 we display an example of 
a functional (in this picture the functional is simply a function) which is not 
differentiable at n . Consequently, there is no unique tangent functional 



at that point. In the example there are in fact infinitely many tangent 
functional at no of which there are two drawn in the plot. 
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Figure 2: The convex functional f 1 has no unique tangent 
functional in the point n . Both 7^ and i 2 are tangent 

functionals 

In the case of infinitely many dimensions one should, however, be careful 
in drawing conclusions from simple plots and we have to resort to formal 
proofs. Since the proof of the previous theorem is important for our discus- 
sion and not difficult to follow we will present it here. Since V is a convex 
set we have that if n e V and m e V then n + A(ni — n ) € V for 
< A < 1 and since F is convex we obtain 

F[n + A(m - n )] < F[n ] + A(F[m] - F[n ]) (210) 

Furthermore since F has a unique continuous tangent functional at no we 
have 

F[n + A(m - n )] > F[n ] + XL[m - n Q ] (211) 



Note that n 1 — n need not be in V but is always in B since that is a linear 
space. If we combine the two inequalities we find that 



FM - F N > F[no + Km-^n )]-F[n ] > ^ _ ^ ^ 

Now since F is convex the function g(X) — F[n + \(rii — n )] — F[n ] is a 
convex real function on the interval [0, 1] and hence continuous. Moreover, 
as the equation above shows, the function g(X)/X has a finite upper and 
lower bound and therefore the limit A — > of g(X)/X exists. We therefore 
find 

F[m] - F[n ] > F'[n , m] > L[m - n ] (213) 

where 

rpu i r F[n Q + X(ni - n )] - F[n ] 

F \n , ni \ — hm — — - (214) 

A^o A 

It remains to show that F'[no,ni] is continuous and linear and equal to L. 
Now F' has the property that for n 2 € V 

rpf, , / mi- F{n + Xe(n 2 - n )} - F{n ] , 

F'[n ,n + e(n 2 - n a )] = hm — ^ ^ ^ ^ = eF'[n ,n 2 ] 

A — *0 A 

(215) 

Therefore from Eq.( 213) we find 

F[n Q + e (m - n )] > F[n ] + eF'[n , m] (216) 

This means that the straight line 

£ = {(n + e(m - no), F[n ] + eF'[n , m])eBx e [0, 1]} (217) 

lies below the graph of functional F. Now since F is convex and lower 
semicontinuous the set of points above the graph of F (the interior points 
of epi(F)) form an open convex set which is nonempty since F is finite 
on V [11]. Now a famous theorem of functional analysis, the Hahn-Banach 
theorem, tells us that if we have an open convex set A (the interior of epi(F) 
in our case) and an affine subspace C (which is the line C in our case) that 
does not intersect A then there is a hyperplane H that contains C in which 
the hyperplane H has the form 

H = {(n,r) e B x K\L'[n] + ar = (3} (218) 

where L' : B — > 1Z is a continuous linear functional and a and (3 are real 
numbers (We can take a ^ since F'[no,ni] is finite and we will have 
no vertical line C or vertical hyperplane 7i). We will not prove the Hahn- 
Banach theorem as it is geometrically intuitive and can be found in most 
textbooks on functional analysis [20, 21]. We will just use its consequences. 



Now we know that our line L is contained in H. The coefficient is then 
determined by the fact that (no, -F[no]) £ H and we find = L'[n ]+aF[no] 
and therefore 

U = {(n, r)eBx n\L'[n - n ] + a{r - F[n }) = 0} (219) 

Now from the fact that (no + e(ni — n ), F[n ] + eF'[no, ni]) eH we find 
that 

L'[n 1 -n o ]+aF'[n o ,n 1 ]= (220) 

From this equation and Eq.(213) we then see immediately that —L'/a is a 
continuous tangent functional at no which then must be equal to L since 
that was the only continuous tangent functional at no- Consequently 

F / [n ,n 1 ] = L[ni-n ] (221) 

and thus 

ii- F[n + \{m - n )] - F[n ] 

— [m - n = lim = L[m - n (222) 

on x^o A 

We therefore obtained what we wanted to prove, F is Gateaux differen- 
tiable at no and SF/5n = L. Now we can apply this theorem to the Lieb 
functional. If we take F — Fl, B = L 1 n L 3 , V — S and use that Fl has a 
unique tangent functional L = — v at every E-V-density and nowhere else, 
we obtain: 

Theorem 10 The functional F^ is Gateaux differentiable for every E- V- 
density in the set S and nowhere else. Moreover the functional derivative 
at an E- V-density is equal to —v where v is the potential that generates this 
density. 

This leaves us with the question which densities in set S actually are 
E-V-densities. A useful result on this point is obtained in the next section. 

14. Ensemble n-representability 

We have seen that the Lieb functional F^ is differentiable at the set of E- 
V-densities in S and nowhere else. For this reason it is desirable to know a 
bit more about these densities. The question therefore is: which densities 
are ensemble v-representable? In this section we will prove a useful result 
which will enable us to put the Kohn-Sham approach on a rigorous basis. 

Theorem 11 The set of E-V-densities is dense in the set S with respect 
to the norm on L 1 n L 3 . 



This statement means the following. Suppose we take an arbitrary 
density no from the set iS, then for every e > we can find an E-V-density 
n such that ||no — n||i < e and ||no — tt.| 1 3 < e. In other words, for every 
density in the set S there is an E-V-density arbitrarily close to it. This 
can also be phrased differently. For every density no m the set S there is a 
series of E-V-dcnsities n k such that 



lim ||n -n fe || p = (223) 

k— >oo 



where the subindex for the norm takes the values p — 1 and p — 3. In order 
to establish this result we need to use a theorem due to Bishop and Phelps. 
For the clarity of the discussion we split the theorem into two parts which 
both yield interesting results for our functional F^. The first part gives 
some insight in the set of potentials V that generate a bound ground state. 
From the relation 

E[v]= inf {F L [n] + / d 3 rn(r)v(r)} (224) 
nes J 

we obtain immediately that for any n E S and any v E L § + L°° 

F L [n] >- J d 3 rn(r)v(r) + E[v] (225) 

Now the functional on the right hand side of the inequality sign is, for 
a given v, a linear functional of n. The inequality sign tells us that this 
functional lies below the graph of F^[n]. A linear functional with this 
property is called F^-boundcd. Let us give a general definition of these 
linear functionals. Let F be a functional F : B — > 1Z from a normcd 
function space (a Banach space) B to the real numbers. Let B* be the 
dual space of B, i.e. the set of continuous linear functionals on B. Then 
L e B* is said to be F-bounded if there is a constant C such that for all 
n e B 

F[n] >L[n] + C (226) 

where the constant C may depend on L but not on n. We therefore sec 
from Eq. (225) that every v e L 2 + L°° defines an F^-boundcd functional 
where L is defined as in Eq.(203). The geometric picture of an F-bounded 
functional is displayed in Figure 3. 




Figure 3: Both L\ and L2 are ^-bounded by the convex 

functional F 



Note that the tangent functionals of Fl are special cases of F^-boundcd 
functionals, they are the Fl bounded functionals that also touch the graph 
of Fl ■ Therefore the set of tangent functionals is a subset of the set of Fl- 
bounded functionals. This can be illustrated with an example. Consider 
the function f(x) = exp (x) on the real axis. This function is convex and 
all tangent functions are of the form g a (x) — ax + (3 where a > and 
(3(a) < 1. The constant function go(x) — (3 with (3 < is an /-bounded 
function, but not a tangent. However, we can always find a tangent with a 
slope arbitrary close to zero, i.e. arbitrary close to the slope of the constant 
/-bounded function go. The following theorem (due to Bishop and Phelps) 
ensures that a similar situation occurs for the case of general Banach spaces: 

Theorem 12 (Bishop-Phelps I) Let F : B — > 1Z be a lower semicontin- 
uous convex functional on a real Banach space B. The functional F can take 
the value +00 but not everywhere. Then the continuous tangent functionals 
to F are B* -norm dense in the set of F -bounded functionals in B* . 



This means that if L is some F-bounded functional then we can find a 
set of tangent functional Lk such that 

lim ||Lo-Lfc|| B . = (227) 

k^oo 

where the limit is taken in the norm on B* . We will not prove this theorem 
here. The proof is clearly described with a geometric interpretation in 
reference [22]. Let us apply this theorem to our functional F^. We know 
that every tangent functional of Fl can be identified with —v, where v is 
a potential that yields a bound ground state, i.e. v G V. Now we know 
that every element — vq for an arbitrary vo G + L°° corresponds to an 
^-bounded functional. Therefore for such a vq there is a series of Vk G V 
such that 

lim Hvo-tfcll = (228) 

k^oo 

where the norm is the L§ +L°°-norm described in Eq.(24). Therefore every 
potential t>o G Li + L°° can be approximated to arbitrary accuracy by a 
potential that yields a bound ground state. This may seem counterintuitive 
at first sight as one can imagine i; to be a repulsive potential, until one 
realizes that one may put the system in big box with a wall of finite height. 
This system will have a bound state if the size of the box is chosen to be 
big enough. The particles will then spread out over the box in order to 
minimize the repulsion between them. If we choose = vq + Wk where Wk 
describes a series of boxes of increasing size but with decreasing height of 
the potential wall then we see that we have created a series of potentials 
in V that approaches wo to an arbitrary accuracy. Let us now discuss the 
second part of the Bishop-Phclps theorem 

Theorem 13 (Bishop-Phelps II) Let F : B — > 1Z be a lower semicon- 
tinuous convex functional on a real Banach space B. The functional F can 
take the value +oo but not everywhere. Suppose uq G B with F[no\ < oo 
and let Lq G B* be an F-bounded functional. Then for every e > there 
exists n f G B and a functional L f G B* such that 

1. \\L t -L \\ B , <e 

2. F[n] > F[n e ] +L e [n~ n e ] for all n. 

3. e||n e - n ||_B < ^[no] - L Q [n ] - mi neB {F[n] - L [n}} 

For the details of the proof we again refer to [22]. The first two points of 
this theorem are equivalent to the previous theorem. They say that any 
.F-bounded functional Lq can be approximated to arbitrary accuracy by a 
tangent functional L e . The inequality in the third point of this theorem 



allows us to make statements about distances between elements in B, which 
in our case will be densities. Note that the right hand side of this inequality 
has the geometric meaning of being the difference of the vertical distance of 
the functional L and F in n and the shortest possible distance between 
Lq and F. Let us now apply this theorem to the functional Fl and show 
that that every density in S is arbitrarily close to an E-V-dcnsity. We first 
need some preliminaries. From Eq.(224) we see that we can write 

F L [n]= sup {E[v]~ f d 3 rn(r)v(r)} (229) 

vehi+L°° 

If the supremum is attained for some v then n is an E-V-density. This 
follows because then there is a density matrix D[n] that yields density n 
(see Theorem 5) such that 

TrD[n]H v = F L [n] + J d 3 rn(r)v(r) = E[v] = MTrDH v (230) 

The density matrix D[n] must therefore be a ground state density matrix. 
If n is not an E-V-density then the supremum is not attained for any v. In 
any case, for every integer k and any density n £ S we can always find 
some v k such that 

E[v k ] - [ d 3 rn (r)v k (r) > sup {E[v\ - j d 3 rn (r)v(r)} - ^ 

= F L [n ]-^ (231) 

where we note that the series v k does not converge to any v if n is not an 
E-V-density. Furthermore, for any n we have 

E[v k ] = mf {F L [n\ + J d 3 rn{v)v k {v)} < F L [n] + J d 3 rn{v)v k {v) (232) 

which in combination with the previous inequality yields 

F L [n] + J d 3 rn{v)v k {v) > E[v k ] > F L [n ] + J d 3 rn a (r)v k (r) - i (233) 

We are now ready to apply the Bishop-Phelps theorem to the Lieb func- 
tional Fl. We take e = 1 , n e S and let — v k correspond to the F-bounded 
functional Lq of the theorem. According to the theorem we can then find 
a tangent functional — w k (the L f of the theorem) such that 



v k - w k \\ < 1 



(234) 



in the norm on L 2 + L°° . From our previous investigations we know that 
the tangent — w k touches the graph of Fl in an E-V-density rife, i.e. 

F L [n] > F L [n k ] - J d 3 rw k (r){n{r) - n k {v)) (235) 

where w k is the potential that generates density n k . The Bishop-Phelps 
theorem then tells us that 

\\n k -n \\ p < F L [n ] + J d 3 rn (r)v k (r) 

- inf{F L [n] + [ d 3 rn(r)v k (r)} (236) 
nes J 

where the normindex has the values p = 1 and p = 3. Note that the 
infimum in this equation can be taken over S rather than L 1 n L 3 since F L 
is defined to be +00 outside S. Then from Eq.(233) we see immediately 
that 

inf {F L [n] + ( ' d 3 rn{v)v k {v)} > F L [n ] + [ d 3 rn {r)v k (r) - \ (237) 
nes J J k 

and therefore 

IK -no||p < ^ (238) 

for p = 1 and p = 3. Since this equation is true for any k we see that any 
density uq G S can be approximated to any accuracy by the E-V-density 
n k in the norms on L 1 f] L 3 . This proves the theorem in the beginning of 
this section. 



15. The Kohn-Sham approach and noninteracting v- 
representability 

We have now come to the discussion of the central equations which form 
the basis of almost any practical application of density functional theory: 
the Kohn-Sham equations. Kohn and Sham [8] introduced an auxiliary 
noninteracting system of particles with the property that it yields the same 
ground state density as the real interacting system. In order to put the 
Kohn-Sham procedure on a rigorous basis we introduce the functional 

T L [n] = inf TrL>f (239) 

We see that this is simply the Lieb functional with the two-particle inter- 
action omitted. All the properties of the functional Fl carry directly over 
to T L - The reason is that all these properties were derived on the basis 



of the variational principle in which we only required that t + W is an 
operator that is bounded from below. This is, however, still true if we omit 
the Coulomb repulsion W. We therefore conclude that Tl is a convex lower 
semicontinuous functional which is differentiable for any density n that is 
ensemble w-representable for the noninteracting system and nowhere else. 
We refer to such densities as noninteracting E-V-densities and denote the 
set of all noninteracting E-V-densities by B . Let us collect all the results 
for Tl in a single theorem: 

Theorem 14 Tl is a convex lower semicontinuous functional with the fol- 
lowing properties: 

1. For any n e S there is a minimizing density matrix D[n] with the 
property Tl[u] = TrD[n]T. 

2. Tl is Gateaux differentiable at the set of noninteracting E-V-densities 
and nowhere else. 

3. The functional derivative at a noninteracting E- V-density n is given 
by: 

where the potential v s generates the density n in a noninteracting 
system. 

From the last point in this theorem we see that if we want to know if 
a given density n from the set S can be obtained as an E- V-density of 
a noninteracting system we may try to calculate the derivative of Tl for 
this density. If it exists then the derivative yields the potential that we 
were looking for. Now Kohn and Sham [8] assumed that for any density of 
an interacting system there is a noninteracting system that has the same 
density as its ground state. We can now ask the question whether the sets 
of interacting and noninteracting E-V-densitics arc equal. This is currently 
not known. However, we can make a number of useful conclusions. First of 
all, if we apply the Bishop-Phelps theorem to the functional Tl we obtain 
the following result: 

Theorem 15 The set of noninteracting E- V-densities is dense in the set 
S with respect to the norm on L 1 n L 3 . 

This means that for any density n € S there is an noninteracting E-V- 
density arbitrarily close. We can in particular choose n to be an interacting 
E- V-density, i.e. n £ B and find a noninteracting E- V-density arbitrarily 
close. We also know from the Bishop-Phelps theorem applied to Fl that the 
set of interacting E-V-densities is dense in S and therefore for any n e S, in 



particular n £ B , we can find an interacting E-V-density arbitrarily close. 
We therefore conclude 

Theorem 16 The set Bo of noninteracting E- V-densities is dense in the 
set B of interacting E- V-densities, and vice versa. 

If we combine this result with previous theorems we obtain the following 
important consequence for the Kohn-Sham scheme: 

Theorem 17 Suppose n G B is an interacting E-V-density. Then for every 
e > there is a noninteracting E- V-density n e € Bo such that 

1. \\n — n e || p < e for p—1 and p = 3. 

2. The density n e is a ground state ensemble density of a noninteracting 
system with potential 

i.e. we can always set up a Kohn-Sham scheme that produces a given in- 
teracting E- V-density to arbitrary accuracy. 

This theorem tells us that in practice we can always set up a Kohn-Sham 
scheme. We also see that if we want to prove that the sets B and Bo of in- 
teracting and noninteracting E- V-densities are equal, then we have to show 
that the potentials v s ^ in this theorem approach some potential v St0 for 
e — ► in some smooth way. This, however, has not been proven until now. 
In numerical calculations (see the discussion at the end of this section) one 
has indeed always succeeded in obtaining a Kohn-Sham potential for given 
interacting E-V-densitics obtained from accurate Configuration Interaction 
(CI) calculations. In these calculations convergence to a given Kohn-Sham 
potential is sometimes difficult to obtain, but seems to happen in a rather 
smooth way. One might therefore expect that the sets of interacting and 
noninteracting E-V-densitics are in fact equal. In order to encourage further 
work in this field we put it here as a conjecture 

Conjecture 1 The sets B and B of interacting and noninteracting E- V- 
densities are equal, i.e. B = Bo- 

Let us describe a couple of cases for which we know this conjecture to be 
true. Consider a system of two particles. Let the interacting system have 
density n. Then we can construct a noninteracting Kohn-Sham system with 
ground state wavefunction 



*jrs(riCTi, r 2 CT 2 ) = ^=^(r 1 )^(r 2 )(a(cr 1 )/3(a 2 ) - (3(a 1 )a(a 2 )) (240) 



where a and (3 are the usual spin functions with a(^) = — 1 and 

ct{—\) — = 0. This wavefunction has a density n(r) = 2|<p(r)| 2 and 
the Kohn-Sham orbital <p satisfies 

(-^V 2 + v s [n](rf) V (r) - e^(r) (241) 

If we now choose ^(r) = yJip{v)/2 then for an interacting E-V-density n 
the Kohn-Sham Hamiltonian with potential 

v s [n](r)= 1 -^^ + e (242) 

has orbital ip as an eigenfunction. Because the density is positive the orbital 
<f = \Jnj2 has no nodes and must be a ground state orbital. We therefore 
have constructed a noninteracting system with ground state n, which is 
even a pure state density. 

Another case where the conjecture is true is for lattice systems. The lattice 
system is obtained by discretizing the many-body Hamiltonian on a grid, 
i.e. we replace the differential equation by a difference equation as one 
might do to solve the Schrodingcr equation numerically. Then any inter- 
acting E-V-density on the grid is also a noninteracting E-V-density. This 
was proven by Chayes, Chayes and Ruskai [23]. This work may be used 
to prove the general conjecture above if one can prove that the continuum 
limit can be taken in some smooth way. 

Let us finally discuss the exchange-correlation functional which is the cen- 
tral object in any application of density functional theory. We define the 
exchange-correlation functional E xc by 

F L [n] =T L [n] + \l d^r ^^f + E xc [n] (243) 

|n -r 2 | 

Since both Tl and Fl are defined on the set S the exchange-correlation 
functional E xc is also defined on that set. Since F L and T L are differentiate 
respectively on the sets B and Bq and nowhere else, the functional E xc is 
differentiate on B fl Bq and nowhere else. The derivative of Eq.(243) on 
that set is given by 

-»C> = --M + /« V jS| + ^ < 2 «> 

If we define the exchange-correlation potential by 



(245) 



then we see that the Kohn-Sham potential can be written as 



v s (r) - v(r) + J d 3 r'0^ + v xc (r) (246) 

We see that on the set B C\ B we have obtained the usual Kohn-Sham 
equations 

f + V + V H [n] + V xc [n]] |$;) = E \^ t ) (247) 

where i = 1 . . . q runs over the q degenerate ground states of the Kohn- 
Sham system. The density must now be calculated from a ground state 
ensemble 

g 

D s [n] = J2»i\®i)(®i\ (248) 
i=i 

of the Kohn-Sham system and is explicitly given as 

<i 

n(r) = Tr£ s [n]n(r) = ^ ^($ 4 |n(r)|^) (249) 

i=l 

These equations can be written more explicitcly in terms of the Kohn-Sham 
orbitals defined to be the eigenstates of the single particle equation 

(- l -V 2 + v s (v)^ t (v) = e m (v) (250) 

It is then easily seen that any Slater determinant 

W) = \<p il ...<Pi„) (251) 

built out of Kohn-Sham orbitals is an eigenstate of the Kohn-Sham system 
with eigenvalue E$ = + . . . + ei N . However, not every eigenstate is 
necessarily a Slater determinant. For instance, if two Slater determinants 
have the same energy eigenvalue, then a linear combination of them also has 
the same eigenvalue. This leads to a subtle point for Kohn-Sham theory of 
degenerate states. In Eq.(249) we require that the density is representable 
by a ground state ensemble of a nonintcracting system. However, we did 
not require that it must be representable by an ensemble of ground state 
Slater determinants. That one should be careful at this point, has been 
shown by Englisch and Englisch [3] and Lieb [1]. They constructed an 
explicit example of a pure state density of a noninteracting system that 
can not be obtained from a single Slater determinant. We will therefore in 
general have 

|$ i ) = 2ay|D j ) (252) 



were \Dj) is a Slater determinant. We further know that every must 
be a ground state for potential v s . This means that all orbitals below the 
highest occupied level must be occupied. If this were not the case then we 
could lower the energy by transferring an electron from the highest level to 
a lower level which would lower the energy and therefore would not be 
a ground state for potential v s . This also means that if two determinants 
\Di) and \Dj) differ, then this difference must be due to different orbitals 
in the highest level. Because both determinants arc cigenfunctions of the 
Kohn-Sham Hamiltonian it is easily seen that these different orbitals must 
have the same Kohn-Sham eigenvalue. This eigenvalue is the Kohn-Sham 
orbital energy of the highest occupied state which we will denote by /i. 
From these considerations we find that the density is of the general form 

n(r)= l^( r )| 2 + E PklM*M(r) (253) 

where (3ki is an Hermitian matrix with eigenvalues between and 1. This 
means that we can diagonalize this matrix with a unitary matrix Uij. We 
then introduce new orbitals for the highest occupied level: 

^(r) = £t/^(r) (254) 

3 

Since these orbitals are linear combinations of degenerate orbitals, they are 
again eigenfunctions of the Kohn-Sham single-particle Hamiltonian for the 
highest occupied level. We can now write the density as 

n(r)= ]T l^(r)| 2 + 5> fc |0 fc (r)| 2 (255) 

with occupation numbers between and 1. We see that for the case of a 
degenerate ground state we have to solve Eq.(250) together with Eq.(255) 
in which also the coefficients rik must be determined. When do we need 
fractional occupation numbers in practice? In practice one sometimes finds 
that when solving the Kohn-Sham equations with occupations equal to 
one that the converged solution has a hole below the highest level. This 
happens for instance for the iron atom within the local density approxima- 
tion [24]. This is an indication that the true Kohn-Sham density could be 
an E-V-density. In practice the lowest energy state is then obtained by a 
procedure called " evaporation of the hole" . In this method one increases 
the occupation of the hole while descreasing the occupation of the highest 
level until both levels are degenerate [24, 25]. The matter has been in- 
vestigated in detail by Schipper et al. [26, 27]. These authors calculated 
by Cl-methods accurate charge densities for a couple of molecules. Sub- 
sequently they constructed the Kohn-Sham potentials that generate these 



densities. A particularly interesting case is the C 2 molecule [26]. When 
one tries to construct a Kohn-Sham potential that produces the ground 
state density of the C 2 molecule using occupation numbers equal to one, 
i.e. for a pure state, one finds that the solution corresponds to a state with 
a hole below the highest occupied level. This is then not a proper Kohn- 
Sham state since it is not the state of lowest energy for the corresponding 
potential, but an excited state. If one then tries to construct an ensem- 
ble solution by the technique of " evaporation of the hole" , one does find a 
proper ground state ensemble for the Kohn-Sham system. From this we can 
conclude that the ground state density of the C 2 molecule is a noninteract- 
ing E-V-density but not a noninteracting PS-V-dcnsity. This means that 
the extension of density functionals to E-V-densities is not only of theoret- 
ical, but also practical interest. As an illustration we display in Figure 4, 
the exchange-correlation potentials of the C 2 molecule for four bonding dis- 
tances. The potentials are plotted along the C — C bond axis as functions of 
the distance z from the bond midpoint. The corresponding bond distances 
R{C—C) are displayed in the inset. The potentials u^f^^ 10 ^ correspond to 
the exchange-correlation potentials of the improper Kohn-Sham solutions 
with a hole. The corresponding state is a pure state consisting of a sin- 
gle Slater determinant. The potentials v^ c are the Kohn-Sham potentials 
corresponding to a ground state ensemble. Both v^-holc and v^ c give 
an accurate representation of the true density of the molecule, although 
the accuracy attained with v^ c is a bit higher. The potential v^ c has the 
features found in many diatomic molecules, i.e. the usual well around the 
atom, the small intershell peaks and a plateau around the bond midpoint. 
The potential v^f^°^ c on the other hand, is heavily distorted as compared 
to v^c- This may be explained by the fact that the wavefunction of the C 2 
molecule has a strong multideterminantial character. With this we mean 
that the simple Hartree-Fock approximation has a relatively small coeffi- 
cient in the Cl-expansion of the wave function. The distortions in v ^c ^ 10 ^ 
therefore appear to be a price we have to pay for producing the density of 
such a strongly multideterminantial state with a single Slater determinant. 
The use of an ensemble formulation of Kohn-Sham theory seems therefore 
especially relevant in cases of strong electron correlations. Other examples 
are provided by near-degeneracy situations such as in avoided crossings [25] 
and the H 2 + H 2 reaction [27] . A recent discussion on Kohn-Sham theory 
for ensembles can be found in Ref. [28] which also gives the beryllium series 
as an explicit example of systems which are ensemble w-representable but 
not pure state w-representable. 




Figure 4: Exchange-correlation potentials producing the accu- 
rate interacting ground state density of the C2 molecule at various 
bonding distances. The potential v^ c is the exchange-correlation 
potential for a ground state ensemble. The potential 

v PS-hole is 

an exchange-correlation potential that reproduces the same den- 
sity for a single Slater determinant with a hole. 



16. The gradient expansion 



Until now we only considered the formal framework of density functional 
theory. However, the theory would be of little use if we would not be able to 
construct good approximate functionals for the exchange-correlation energy 
and exchange-correlation functional. Historically the first approximation 
for the exchange-correlation functional to be used was the local density 
approximation. In this approximation the exchange-correlation functional 
is taken to be 

E™ A [n] = / d 3 re xc (n(v)) (256) 

where e xc (n) is the exchange-correlation energy per volume unit for a homo- 
geneous electron gas with density n. We therefore treat the inhomogcncous 
system locally as an electron gas. This simple and crude approximation 
turns out to be surprisingly succesful for realistic and very inhomogcncous 
systems although one would expect that the LDA would only be accurate 
for systems with slowly varying densities. The LDA therefore seems a suit- 
able starting point for more accurate approximations. In this section we 
will show that the LDA is in fact the first term in a systematic expan- 
sion of the exchange-correlation functional in terms of spatial derivatives 
with respect to the density. This expansion is commonly referred to as the 
gradient expansion [5, 29, 30] and has the following form 

E xc [n] = E^ A [n] + J rf 3 r Sl (n(r))(V«(r)) 2 

+ y>r 52 (n(r))(V 2 n(r)) 2 + ... (257) 

where the functions gi(n) are uniquely determined by the density response 
functions (of arbitrary order) of the homogeneous electron gas. The gradi- 
ent expansion presents an, in principle, exact way to construct the exchange- 
correlation functional for solids, provided the series converges. We will come 
back to the question of convergence. To derive the gradient expansion we 
start out by expanding the exact exchange-correlation energy functional 
around its homogeneous electron gas value n as follows 

E xc [n] = E xc [n }+ 

oo . 

-r / d 3m rK x ^(n ;v 1 ...v m )Sn(v 1 )...Sn(v m ) (258) 

m=l 

where d 3 r = d 3 r\ . . . d 3 r m and n(r) = n + Sn(r). We further defined 

X m TP 

^ fl - r "' = ^)...M, B ) ^ (259) 



The first term E xc [n ] in Eq.(258) is the exchange-correlation energy of a 
homogeneous system with constant density no and is therefore a function of 
no rather than a functional. We will therefore write this term as E xc (no). 
This function is by now well-known from extensive investigations of the 
homogeneous electron gas [31]. Since the electron gas has translational, 
rotational and inversion symmetry the functions functions K x ™^ satisfy 

tf<?>(no; r!... r m ) = (n ;r 1 +a... r m + a) (260) 

^(no;ri...r m ) - K x f\n ;Rr 1 ...Rr m ) (261) 
(n<>; r!... r m ) = K x ™\n ; -n . . . - r m ) (262) 

where a is an arbitrary translation vector and where i? is an arbitrary 
rotation matrix. Furthermore, the function K x ™^ has full permutational 
symmetry, i.e. 

K x ™\n ; r 1 ...r m ) = K x ™\n ; r P(1) . . . r P(m) ) (263) 

for an arbitrary permutation P of the indices 1 . . . m. If we choose a = — v\ 
in Eq. (260) we see that 

K x ™\n Q ;v 1 ...v m ) = K x ™\n ;0,r 2 - n . . .r m - n) 

= L( m )(n ;r 2 -r 1 ...r m -r 1 ) (264) 

where the latter equation defines the function as a function of m — 1 
variables. Let us give some specific examples. The function K X X } is given 

by 

xW(n ;ri) = fa ^ | n=no = u xc (n ) (265) 

and we see that this function is simply the exchange-correlation potential 
of the electron gas which, due to translational invariance, does not depend 
on ri. Since Sn(r) integrates to zero the term with K^J in the expansion 
Eq.(258) does not contribute to E xc [n\. The first nontrivial term is therefore 

K x 2 J(n ;r u r 2 ) = __-^__| n=Ilo = L^(n ; r 2 - n) (266) 

which leads to an expansion of the form 

E xc [n} = E xc {n ) + ^ J d 3 r 1 d 3 r 2 L^(n ;r 2 -r 1 )Sn(r 1 )8n(r 2 ) + . . . (267) 

The function (often denoted by K xc in the literature) has been the 
subject of many investigations [32, 33]. Let us now go back to the more 
general case. If we introduce the Fourier transform of a function / as follows 

/(qi . . . q m ) = / rf 3m r/(n . . . r m ) e -*ii-'i---*im-'m (26 8) 



then 

&&°(no;qi . . . q m ) = (27r) 3 <5( qi + . . . + a m )U m \n a ; q 2 . . . q m ) (269) 
The expansion for E xc then becomes 
E xc [n] = E xc (n )+ 

°° 1 f d 3m CI 

53 ^ / (2^*~ ^ * ■ • ■ 1™)*»(«U) • • ■ (270) 
m— 1 ^ ' 

where d 3m q = d 3 qi . . .d 3 q m . This can also be written in terms of the 
functions as follows 



E xc [n] = E xc (n )+ 

EJ_ f d 3 q 2 d 3 q,, 
ml J (2tt) 3 "' (27r) 3J 

x<5n(-q 2 - ... - q m )<5n(q 2 ) . . . Sh(q m ) (271) 



r L( m )(n ;q 2 ...q m ) 



The gradient expansion is then obtained by expanding in powers of 

around — and subsequently transforming back to real space. If we 
do this our final expansion will then still depend on our reference density 
n which is undesirable for application of the functional to general systems. 
However, we will show that we can get rid of this dependence on no by 
doing infinite resummations. The key formula that will enable us to do 
these resummations is the following first order change of K x ^ : 

(ri . . .r m ) = J tPrK£ +1 \no;T, ri . ..r ro )Jn(r) (272) 

If we take Sn(r) — Sn to be a constant change we obtain 

dK x m) 



dn 



-(no;ri...r m ) = J d 3 rK t +1 \n - r, n . . . r m ) (273) 



Note that by taking <5n(r) = Sno we violate the condition that 6n(r) inte- 
grates to zero. Nevertheless Eq.(273) can alternatively be derived directly 
from the properties of the response functions [30] and is related to a gener- 
alized form of the well-known compressibility sumrule of the electron gas. 
An important special case is obtained for m = 1: 

d 3 r 2 L^(n ; r 2 - n) = L^(q = 0) (274) 



We see that the q = value of can be directly calculated from the 
knowledge of v xc (no). For m > 2 relation (273) directly implies the follow- 
ing relation between L^> and L^ m+1 ' 

— (n ;q2 •■•q m ) = L (m+1) (n ; -q 2 - . . . - q m , q 2 . . . q m ) (275) 

Cno 

The importance of this formula will become clear if we work out, as an 
explicit example, the lowest order in the gradient expansion. From the 
symmetry properties of KxT^ one finds that its Fourier transform has the 
following expansion in powers of q 4 : 

KM (no; qi . . . q„) = (2 7 r) 3 ( 5(q 1 + . . . + q m ) 

x (Kt ] + K[ m) p[ m \ qi . . . q m ) + K^PtX^ . . . q m ) + . . .) (276) 

in which the K^™^ are coefficients and are symmetric polynomials. 

This is a direct consequence of the permutational symmetry of K xc . The 
first polynomial p[ m ^ (taking into account the 5-function) has the explicit 
form 

P 1 (m) ( qi ... qm ) = q ? + ... + q ^ (277) 

and is invariant under permutation of the indices, i.e. it transforms accord- 
ing to a one-dimensional representation of the permutation group. More 
details on the group theoretical treatment of these reponse functions can be 
found in the work of Svendsen and von Barth [34] . From the properties of 
we find that for m > 2 the function has the following expansion 

i (m) (no; q 2 • • • q m ) - L^ n \n ) + '(n )P (m) '(<to • • • «k) + • • ■ (278) 
where the polynomial is given by 

P (m) (q2---q m ) = (-q 2 -...-q m ) 2 + q2 + --- + q™ 

m m 

= 2^q, 2 -2^ q;. qj (279) 

i=2 i>j>2 

It is easily seen that this polynomial has the property 

p(™+i)(-q 2 - . . . - q m , qa . . . q m ) = p^(q 2 • • • q m ) (280) 

If we use this property together with Eq.(278) and Eq.(275) we obtain for 
the coefficients in the expansion of L^" 1 ) the following equations 

ft T ( m ) 

= L o m+1) ("°) ( 281 ) 

= L ( r +1 \n ) (282) 



where to > 2. Together with 

* m <* = <» = feM = <*»> 

where e xc (no) is the exchange energy per volume unit of the electron gas, 
this yields 

4 m) M = ^p(no) (284) 

^Vo) = ^prK) ( 2 85) 

It is these two equations that will allow us to eliminate the dependence on 
no in the lowest order gradient expansion. If we insert the explicit form of 
f Eq.(278) into Eq.(271) and Fourier transform back to real space we 
obtain 

E xc [n] = E xc (n a )+ 

oo - 

E^i L o m) («o) d\{8n{v)r + 

oo 

V — L[ m \n Q )m(m - 1) / d 3 r(VSn(r)) 2 (Sn(r)) m - 2 + . . .(286) 



m=2 



Since V<5n(r) = Vn(r) this can be rewritten with help of Eq.(284) and 
Eq.(285) as 



E xc [n] = E xc (n )+ 

y 1 f rf 3 r(Mr)r + 

^ to! dnV? J K y " 



m=l """ 

-2r(2) 



E (m l 2) , g^J ° } / rfMVn(r))^(Mr)) m - 2 + ...(287) 

We see that the first two terms simply give the expansion of the LDA 
functional around the constant density n ,i.e. 

E x ? A [n] = J d 3 re xc (n + Sn(r)) (288) 

whereas the third term involves an expansion of the coefficient L^ 2 \n a + 
5n(r)) around uq. We obtain 

E xc [n] = E™ A [n] + J d 3 rL {2) (n(r))(Vn(r)) 2 + . . . (289) 



and wc sec that we succeeded in eliminating the dependence on n in the 
gradient expansion. Morever we see that the local density approximation 
naturally appears as the first term in the gradient expansion. Furthermore 
we find that the coefficient gi(n) in our first equation (257) is completely 
determined by the q-expansion of the function L^(no; q) to order q 2 . We 
also see that the replacement n — > n(r) requires a resummation over re- 
sponse functions of arbitrary order. The dependence on n can also be 
removed to higher order in powers of q. For an explicit example up to 
order q 4 we refer again to Svcndsen and von Barth [34]. 
Let us now address the question of convergence of the gradient expan- 
sion. This question has been investigated for the gradient expansion of the 
exchange-energy functional for which a comparison with exact exchange 
energies is possible. For this case the analytic form of ZA 2 )(no;q) is known 
from the impressive work of Engel and Vosko [35, 36] and some higher order 
gradient coeficients have been determined by Svendsen, Springer and von 
Barth [34, 37]. One finds [37] that the gradient expansion performs very 
well for metallic systems, but that the quality deteriorates as soon as the 
system acquires an energy gap. This may not be so surprising if one real- 
izes that the appearance of a gap drastically changes the low q-behavior of 
the response functions which determine the gradient coefficients [38]. 

This means that the standard gradient expansion can not deal with insu- 
lators or finite systems (which may be modelled as insulators by repeating 
them periodically with a large lattice constant). Further progress along the 
lines of the gradient expansion may be obtained by study of the so-called 
gapped electron gas [39]. Currently the most fruitful approach to the con- 
struction of simple and accurate functionals is the so-called Generalized 
Gradient Approximation (GGA) [40, 41]. In this approach one specifies a 
form of the pair-correlation function of the many-electron system and de- 
termines the parameters in this function by sumrules and information from 
the straightforward gradient expansion. These functionals have lead to 
large improvements in molecular binding energies as compared to the local 
density approximation [42]. However, since the approach is not systematic 
it is difficult to improve the quality of the current GGA functionals. 

17. The optimized potential method and the e 2 -expansion 

In this section we will describe a second systematic method to construct 
density functionals, namely perturbation theory starting from the Kohn- 
Sham Hamiltonian [43, 44]. The method is based on traditional pertur- 
bation theory and is comparable in computational cost. For this reason 
the method is less suited to the calculation of properties of large systems. 
However, the method has the theoretical advantage that it can be used as a 



benchmark to test the quality of different approximate density functional. 
We consider the Hamiltonian 



H e = H s + e(H-H s ) = H s + e(W-V Hxc ) (290) 

where 

* = ^ (291) 
i>j |r * r ^ 

AT 

Vhxc = Y^VHfa) + Vxc (Ti) (292) 

where e 2 is the square of the electronic charge and vh is the usual Hartree 
potential 

„2 f ,3 / ™( r ') 



Mr) = e ^dV^^ (293) 

and v xc the exchange-correlation potential. The perturbation term is simply 
the difference between the true and the Kohn-Sham Hamiltonian and we 
are interested in the case e = 1 although this makes expansion in powers 
of e rather doubtful. We will come back to the point of convergence later. 
We can now do standard perturbation theory and obtain the ground state 
energy 

E(e) = E S + eE {1) + e 2 E {2) + ... (294) 

where E s is ground state energy of the auxiliary Kohn-Sham system, i.e. 
simply a sum over orbital energies and the terms E^ is the energy to order 
i in powers of e. The first two terms are explicitcly given by [43, 44] 

= (<S>s\W-V Hxc \$> s ) (295) 

Em = g l(^-vwi<MI 2 (296) 

i— 1 s s,z 

where E s ^ and <fr s ^ are the excited state energies and wave functions of 

the Kohn-Sham system. We now note that the energies are implicit 

functionals of the density through their dependence on the Kohn-Sham 
potential, orbitals, and energies, i.e. 

E^[n] = EM[{y k [n]},{e k [n}},v xc [n]\ (297) 

This follows directly from the Hohenberg-Kohn theorem applied to a non- 
interacting system. The density n uniquely determines the Kohn-Sham 
potential v s (up to a constant) and therefore the also the orbitals (up to 
a phase factor) and eigenvalues (up to constant). The arbitrariness with 



respect to a constant shift and with respect to the phase factor cancels out 
in the energy expression and therefore the i-th order energy becomes a pure 
density functional. We therefore have the following series of implications 

n(r) - v,(t) - {p fc (r), e fe } - E« (298) 

Note that the perturbation theory that we constructed is not yet in a form 
that we can use in practical calculations. This is because the perturbing 
Hamiltonian contains the exchange-correlation potential which is unknown 
from the start and has to be determined sclf-consistcntly from the pertur- 
bation series. However, the equations can be simplified if we expand the 
energies and potentials in powers of the interaction strength e 2 and take 
e=l. This leads to the following set of equations [44] 

E[n] = T s [n] +j d 3 rn(r)v(r) 

+ y/^'t^ + E^ ( 2 ") 

' ' 2 = 1 

OO 

v Hxc (r) = v H (r) + J2e 2i v£(r) (300) 
i=i 

45M = jgf (3oi) 

The e 2 and e 4 terms in the expansion of the exchange-correlation energy 
have the explicit form 

e 2 £«N = <* a |W|*.) -^J dW^^l (302) 



and 



[lBa ^Mz|z|fM ,303, 

i=l 

The term of order e 2 can be written explicitly in terms of the Kohn-Sham 
orbitals as follows: 

4'»H —\tf W/ ^ff"" (304) 

We see that this expression has exactly the same form as the usual expres- 
sion of the exchange energy within the Hartrcc-Fock approximation. The 
difference, however, is that the orbitals in this expression are Kohn-Sham 



orbitals which, in contrast to the Hartree-Fock orbitals, are eigenfunctions 
of a single-particle Hamiltonian with a local potential. The numerical value 
of K^jn] and the Hartree-Fock exchange will therefore in general differ 
from each other. However, because of the similarity to the Hartree-Fock 
definition of exchange we will define the exchange functional within density 
functional theory to be E x [n] = E^Jln]. Corresponding to the exchange- 
functional there is a local exchange potential defined as 

We see from Eq.(303) that we need to know this potential in order to 
calculate the e 4 contribution to the exchange-correlation energy. This is 
a general feature of the present perturbation theory. In order to calculate 
Exc we need to calculate vie ^ first. Let us therefore start by calculating 
v x . For this we have to calculate the first order change SE X in the exchange 
functional due to a change Sn in the density. Since densities and potentials 
are in 1-1-correspondence this task amounts to calculating the change in 
E x due to a change Sv s in the Kohn-Sham potential. This is readily done 
by perturbation theory. The change S<p k and Se k of the orbitals and orbital 
energies of the Kohn-Sham system due to a small change 5v s in the potential 
is given by solution of the equation 



1 V 2 + v s (r) - e k I ^ fe (r) = (Se k - Sv s (r))<p k (r) (306) 



2 



This gives 



Se k = J d 3 r<pi(r)cp k (r)5v s (r) (307) 

tyfc(r) = - j d?r'G k {v,v')^ k {v')Sv s {v') (308) 

GhM = E ^MM (309) 

From these equations we obtain the following functional derivatives 

' )(k = ^W^(r) (310) 

= -G fc (r,i> fe (r') (311) 



8v s (r) 
S^ k (r) 



Sv s (r>) 

The density change Sn is given by 

N 



Sn(r) = vlWMr) + c - c - = I d 3 r' Xs (r, r')Sv s (r') (312) 
fe=i J 



where we defined the static density response function x.s of the Kohn-Sham 
system by 



Xs(r ' r ' } = J^) = ~ 2 vl(r)G k (r, r')^(r') + c.c. 



(313) 



k=i 



We can now readily derive the following integral equation for v x : 

5E X 



A x (r) = 



5u s (r) 



d 3 r / fa(rQ 
Sn(r') 5v s (r) 



= J dV Xs (r,r> x (r') (314) 



where the inhomogeneity is given by 



N r t)F 

Ax(r) = -^y dV^-^G fc (r',r)^(r)+ C .c. (315) 

Since both and are given as explicit functionals of the Kohn-Sham 
orbitals and orbital energies, the exchange potential can be found from 
simultaneous solution of the equations 



efe^fcW 



1 



V z + v(r) + v H (r) + eV(r) <p k ( T ) ( 316 ) 



A,(r) = J dV Xs (r,r'K(r') 



(317) 



Once we have solved these equations and determined v x we can go on 
and calculate E X 2 J from Eq.(303). To calculate the e 6 contribution to the 

(2) 

exchange-correlation energy we first have to evaluate v xc . This potential 
is the solution of the integral equation 



(318) 



(2) 

Since E xc is an explicit functional of the orbitals, orbital energies and v x 
the inhomogeneity K x 2 } can be calculated from 

Ai 2 2(r) = 



+ 



k=l 



CO „ 

I — 1 J 



d A r 



k=i 



6ip k (r') Sv s (r) 

j 5E X 2 J 5e k 
8e k 5v s ( 



d A r 



, 5E X 2 } 5v x {v') 
5v x (v') Sv s (r) 



(319) 



All terms in this equation are explicitly known, except for the term 5v x /Sv s . 
However, for this term we can find an integral equation by differentiation 



of Eq.(317) with respect to v. 



^Ax(ri) f , 3 <5x s (ri,r 3 f 3 ^x(r 3 ) , . 

x / x = / rf r 3 x 7 x ^x(r 3 )+ / d J T- 3 x s (ri,r 3 ) (320 
0v s (r 2 ) J ov s (r 2 ) J Sv s (r 2 ) 

Since both A x and Xs are explicitly known in terms of orbitals and orbital 
energies their derivatives with respect to v s are also explicitly known in 
terms of these quantities (sec reference [44] for more details) and therefore 
Eq.(320) determines Sv x /Sv s uniquely So we see that in order to deter- 

(2) 

mine v xc we have to solve two integral equations. For realistic systems 
these equations have only be solved approximately (for an explicit solution 
of Eq.(320) see reference [34])). From v xc one could go on along the same 
lines to determine E xc ' and subsequently v xc . The determination of these 
higher order energies and potentials becomes more involved. Nevertheless, 
the perturbation scries represents an explicit construction procedure for the 
exact exchange-correlation energy and potential, provided that the series 
converges. Before we go on to discuss the convergence properties of this 
scries, let us briefly review the first order equations from a different view- 
point. It is readily seen that one can write the total energy to order e 2 , 
which we denote by £a[n], as 

E^n] = T s [n] + J d\n{v)v{v) + & - J d 3 rd 3 r' n ^ H ^ + e 2 E x [n] 

= ($s\H\$ s ) (321) 

This is just the expectation value of the true Hamiltonian H of the system 
with the Kohn-Sham wave function |$ s ). Because of the 1-1-correspondcncc 
between the density n and the potential v s we can also regard E\ as a 
functional of the potential v s , i.e. 

Ei[v.] = {* a [v a ]\H\* a [v a ]) (322) 

We may now try to find an approximation to the true total energy of the 
system by choosing a local potential v s that minimizes the energy expression 
Ei [v s ] . This means that we have to solve the variational equation 



°-6v.(T)-£>J dr 6<p k (r>) 5v s (r) +CX - 



fe=i 



If one works out this expression one obtains equations that are identical 
to Eqns.(316) and (317). These equations were first derived by Talman 
and Shadwick [45]. Since in our procedure we optimized the energy of 
a Slater determinant wavefunction under the constraint that the orbitals 



in the Slater determinant come from a local potential, the method is also 
known as the Optimized Potential Method (OPM). We have therefore ob- 
tained the result that the OPM and and the expansion to order e 2 are 
equivalent procedures. The OPM has many similarities to the Hartree- 
Fock approach. Within the Hartree-Fock approximation one minimizes the 
energy of a Slater determinant wavcfunction under the constraint that the 
orbitals are orthonormal. One then obtains one-particle equations for the 
orbitals that contain a nonlocal potential. Within the OPM, on the other 
hand, one adds the additional requirement that the orbitals must satisfy 
single-particle equations with a local potential. Due to this constraint the 
OPM total energy E\ will in general be higher than the Hartree-Fock en- 
ergy Ejfp, i.e. Ei > Eh p. We refer to [46, 47] for an application of the 
OPM method for molecules. 

We finally make some comments on the calculation of functional derivatives 
in this section. We stress this point since careless use of the chain rule for 
differentiation has led to wrong results in the literature [48] . As an example 
we consider the exchange functional E x . When we regard this functional 
as an explicit functional of the orbitals then the functional derivative with 
respect to ipk is given by 

(a>)-g/ A 'f^'> <™ 

where we used the subindcx e to indicate that we regard the functional 
as an explicit orbital functional. This functional derivative represents the 
change in E x due to a change Stpk in orbital ifk, while keeping all other 
orbitals fixed. Moreover, we regard ip^ and ip* k as independent variables. If 
we regard E x as a density functional we must require that all orbitals arc 
cigenfunctions of a nonintcracting Hamiltonian with a local potential v s . It 
is clear that we can we can never find a change 5v s in a local potential that 
induces a change in only one orbital while keeping the other orbitals fixed. 
If a potential changes one orbital tpk to tpk + Stfk then all other orbitals will 
change too. The change in E x regarded as a functional of v s is then given 
by the functional derivative 

SE X A r 3 , / SE X \ 5ip k (r') 




(325) 



We see that in this case the explicit orbital derivatives only occur in the 
sum in which all the orbital changes must be taken into account. We may 



also introduce an implicit derivative (SE X /Sipk)i, which can be given the 
meaning of giving the change in E x if we know that there is a potential 
change Sv s that changes ifk to ipk + Sipk- This implicit orbital derivative, 
which we denote by subidex i is expressed in terms of the explicit orbital 
derivatives as follows 

SE X \ _ y, r 3 ^, ( SE X \ 5ipi(r') 



i=l 
N 

E 



6( Pk(r)J l f-'J \5<Pi{r') ) 5tp h (T) 



( ^] ^ (326) 
* 6<p*(r>)J e ^ fe (r) 



The derivatives 5ipi/5ipk appear now to take into account that, if orbital (fk 
changes to ipk + fitfik by some potential change Sv s , then the other orbitals 
ipi change to (fi + Sipi. By using the implicit orbital derivative we regard 
E x as a functional of the potential, or equivalently of the density, and we 
can therefore use the chain rule 



5E. 



In this equation (Sn/Sifk)i is also an implicit derivative given by: 



6n(r') 



N 



.^ (r) ^W + ^ (r) ^) 



(328) 



which should be compared to the explicit derivative 

' <5n(r') 



<fy>k(r) 



= <P%(t) (329) 



for k = 1 . . . N. We stress that Eq.(327) is only true when using the im- 
plicit derivatives. If in Eq.(327) we replace the implicit derivatives by the 
explicit ones of Eq.(324) and Eq.(329) we obtain a result that is wrong. As 
pointed out this is due to the fact that a change in just one orbital can 
not be induced by a change in a local potential which in density functional 
theory is in 1-1-correspondence with the density. Most of the confusion can 
be avoided by regarding all functionals as functionals of the potential v s 
and by calculating the change in the functionals by means of perturbation 
theory. In this way one can avoid use of the implicit orbital derivatives as 
in Eq.(326). 

Let us finally come back to the question of convergence of the perturbation 
series. The perturbation series presented in this section is very similar to 
M0ller-Plesset perturbation theory starting from the Hartree-Fock approx- 
imation. For the M0ller-Plesset perturbation theory it is known that it is 



in general divergent [49, 50, 51]. However, it is well-known that when car- 
ried out to low orders this perturbation theory gives reasonable answers. 
The M0ller-Plesset perturbation series has therefore all the features of an 
asymptotic series [52]. Since the perturbation series in this section is very 
similar to the M0ller-Plesset scries it will in general also diverge. This 
has indeed been found in a perturbation theory on the basis of some ap- 
proximate Kohn-Sham Hamiltonians [53]. Nevertheless, it has been found 
for the method presented in this section that low orders in perturbation 
theory give good results [54] and we conclude that our series give at least 
an asymptotic expansion for the exchange-correlation energy and potential. 



18. Outlook and conclusions 

In this work we have given on overview of the mathematical foundations 
of stationary density functional theory. We discussed in great detail the 
question of differentiability of the functionals and showed that the Kohn- 
Sham theory can be put on a solid basis for all practical purposes, since 
the set of noninteracting E-V-densities is dense in the set of interacting 
E-V-densities. The question whether these two sets are in fact identical is 
still an open question. We further discussed two systematic approaches for 
the construction of the exchange-correlation functional and potential. 
What can we say about future developments within density functional the- 
ory? There have been many extensions of density functional theory in- 
volving spins, relativistic effects, temperature, superconductivity and time- 
dependent phenomena. The last few years we have, for instance, seen many 
applications of response properties, rather than ground state properties, us- 
ing time-dependent density functional theory. For these extended density 
functional theories it is, of course, more difficult to provide a rigorous the- 
oretical basis. This is, however, not particular to density functional theory, 
but applies to any method that deals with many-body systems, especially if 
one is interested in phenomena like superconductivity or interactions with 
laser fields. Nevertheless, also for these complicated cases simple density 
functionals have provided encouraging results although there is still a clear 
need for more accurate density functionals. In view of the success of den- 
sity functional theory for ground state calculations it seems worthwhile to 
explore these new areas. 
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